Рассчитать площадь под кривой

Я хотел бы рассчитать площадь под кривой, чтобы сделать интеграцию без определения функции, такой как в integrate(),

Мои данные выглядят так:

Date          Strike     Volatility
2003-01-01    20         0.2
2003-01-01    30         0.3
2003-01-01    40         0.4
etc.

Я построил plot(strike, volatility) смотреть на изменчивую улыбку. Есть ли способ интегрировать эту построенную "кривую"?

Источник

7 ответов

Решение

AUC аппроксимируется довольно легко, если посмотреть на множество фигур трапеции, каждый раз между x_i, x_{i+1}, y{i+1} а также y_i, Используя rollmean пакета zoo, вы можете сделать:

library(zoo)

x <- 1:10
y <- 3*x+25
id <- order(x)

AUC <- sum(diff(x[id])*rollmean(y[id],2))

Убедитесь, что вы заказываете значения х, иначе ваш результат не будет иметь смысла. Если у вас есть отрицательные значения где-то вдоль оси y, вам нужно выяснить, как именно вы хотите определить область под кривой, и отрегулировать соответственно (например, используя abs())

Что касается ваших последующих действий: если у вас нет формальной функции, как бы вы ее составили? Так что, если у вас есть только значения, единственное, что вы можете приблизить, это определенный интеграл. Даже если у вас есть функция в R, вы можете вычислить только определенные интегралы, используя integrate(), Построение формальной функции возможно только в том случае, если вы также можете определить ее.

Источник

Просто добавьте следующее в вашу программу, и вы получите площадь под кривой:

require(pracma)
AUC = trapz(strike,volatility)

От ?trapz:

Этот подход в точности соответствует приближению для интегрирования функции с использованием правила трапеции с базовыми точками x.

Источник

Еще три варианта, в том числе один с использованием метода сплайна и один с использованием правила Симпсона...

# get data
n <- 100
mean <- 50
sd <- 50

x <- seq(20, 80, length=n)
y <- dnorm(x, mean, sd) *100

# using sintegral in Bolstad2
require(Bolstad2)
sintegral(x,y)$int

# using auc in MESS
require(MESS)
auc(x,y, type = 'spline')

# using integrate.xy in sfsmisc
require(sfsmisc)
integrate.xy(x,y)

Трапециевидный метод менее точен, чем сплайн-метод, поэтому MESS::auc (использует сплайн метод) или Bolstad2::sintegral (использует правило Симпсона) должно быть предпочтительнее. Их DIY-версии (и дополнительный подход с использованием квадратурного правила) находятся здесь: http://www.r-bloggers.com/one-dimensional-integrals/

Источник

Хорошо, так что я приезжаю немного поздно на вечеринку, но просматриваю ответы просто R Решение проблемы отсутствует. Здесь все просто и чисто:

sum(diff(x) * (head(y,-1)+tail(y,-1)))/2

Решение для OP тогда читается как:

sum(diff(strike) * (head(volatility,-1)+tail(volatility,-1)))/2

Это эффективно вычисляет площадь, используя трапециевидный метод, принимая среднее значение "левого" и "правого" значений y.

NB: как уже отметил @Joris, вы можете использовать abs(y) если это будет иметь больше смысла.

Источник

В мире фармакокинетики (ПК) вычисление различных типов AUC является обычной и фундаментальной задачей. Существует множество различных расчетов AUC для фармакокинетики, таких как

  • AUC0-t = AUC от нуля до времени t
  • AUC0-last = AUC от нуля до последнего момента времени (может быть таким же, как указано выше)
  • AUC0-inf = AUC от нуля до бесконечности времени
  • AUCint = AUC за промежуток времени
  • AUCall = AUC за весь период времени, за который существуют данные

Один из лучших пакетов, который делает эти вычисления, является относительно новым пакетом PKNCA от людей в Pfizer. Проверьте это.

Источник

Ответ Joris Meys был великолепен, но я изо всех сил пытался удалить NA из моих образцов. Вот небольшая функция, которую я написал, чтобы справиться с ними:

library(zoo) #for the rollmean function

######
#' Calculate the Area Under Curve of y~x
#'
#'@param y Your y values (measures ?)
#'@param x Your x values (time ?)
#'@param start : The first x value 
#'@param stop : The last x value
#'@param na.stop : returns NA if one value is NA
#'@param ex.na.stop : returns NA if the first or the last value is NA
#'
#'@examples 
#'myX = 1:5
#'myY = c(17, 25, NA, 35, 56)
#'auc(myY, myX)
#'auc(myY, myX, na.stop=TRUE)
#'myY = c(17, 25, 28, 35, NA)
#'auc(myY, myX, ex.na.stop=FALSE)
auc = function(y, x, start=first(x), stop=last(x), na.stop=FALSE, ex.na.stop=TRUE){
  if(all(is.na(y))) return(NA)
  bounds = which(x==start):which(x==stop)
  x=x[bounds]
  y=y[bounds]
  r = which(is.na(y))
  if(length(r)>0){
    if(na.stop==TRUE) return(NA)
    if(ex.na.stop==TRUE & (is.na(first(y)) | is.na(last(y)))) return(NA)
    if(is.na(last(y))) warning("Last value is NA, so this AUC is bad and you should feel bad", call. = FALSE) 
    if(is.na(first(y))) warning("First value is NA, so this AUC is bad and you should feel bad", call. = FALSE) 
    x = x[-r]
    y = y[-r]
  }
  sum(diff(x[order(x)])*rollmean(y[order(x)],2))
}

Затем я использую его с приложением на моем фрейме данных: myDF$auc = apply(myDF, MARGIN=1, FUN=auc, x=c(0,5,10,15,20))

Надеюсь, что это может помочь новичкам, как я:-)

РЕДАКТИРОВАТЬ: добавлены границы

Источник

Вы можете использовать пакет ROCR, где следующие строки дадут вам AUC:

pred <- prediction(classifier.labels, actual.labs)
attributes(performance(pred, 'auc'))$y.values[[1]]
Источник
Другие вопросы по тегам