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

Question

Я хотел бы рассчитать площадь под кривой, чтобы выполнить интеграцию без определения функции, например, в 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), чтобы посмотреть на улыбку волатильности. Есть ли способ интегрировать эту построенную «кривую»?

Answer 1

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))

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

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

Answer 2

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

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

От ?trapz:

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

Answer 3

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

# 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 (использует правило Симпсона), вероятно, предпочтительнее.Сделай сам их версии (и дополнительный подход, использующий правило квадратуры) здесь: http://www.r -bloggers.com / одномерные интегралы /

Answer 4

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

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

Решение для OP затем читается как:

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

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

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

Answer 5

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

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

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

Answer 6

Ответ 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))

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

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

Answer 7

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

pred <- prediction(classifier.labels, actual.labs)
attributes(performance(pred, 'auc'))$y.values[[1]]