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

Ответ 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(). Построение формальной функции возможно только в том случае, если вы также можете определить ее.

Ответ 2

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

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

От ?trapz:

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

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

Ответ 4

ОК, поэтому я немного опаздываю на вечеринку, но, пройдя ответы, простое решение 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), если это имеет смысл.

Ответ 5

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

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

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

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

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

EDIT: добавлены ограничения

Ответ 7

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

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