logo标签列表

计算曲线下的面积

rnumerical-integration
71

我想计算曲线下的面积以进行积分，而不需要定义一个像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)以查看波动率笑曲线。有没有一种方法可以对这个绘制的“曲线”进行积分？

1
请看这个相关的问题：https://dev59.com/4m445IYBdhLWcg3wWZAU - Andrie
3
@Andrie: 那是一种不同类型的AUC... - Joris Meys
7个回答
55

通过查看许多梯形图形的面积，可以很容易地近似计算AUC，每次限于x_i，x_{i+1}，y{i+1}和y_i之间。 使用zoo包的rollmean函数，您可以执行以下操作：

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()计算定积分。仅当您还可以定义该函数时，才可能绘制正式函数。

这很棒，但如果有一些值缺失，公式就不再起作用了。 - Dan Chaltiel
@DanChaltiel 如果有一些值缺失，就无法知道曲线下面积的真实值。所以这对我来说似乎不是问题。如果你想忽略缺失的数据，在计算之前只需删除缺失的观测值即可。 - Joris Meys
如果你有10个x值，但只有9个y值，如果你不计算缺失值，你可以相当准确地近似AUC。对于我来说，删除所有只有一个NA的样本似乎是一种浪费。 - Dan Chaltiel
这确实是个问题，因为作为R的新手，我不得不花费一些时间才能找到如何删除NA值（我刚刚在这里发布了它）。请不要忘记R并不是很容易学习的。你觉得容易的东西对于每个人来说都不会容易。 - Dan Chaltiel
@DanChaltiel，这就是为什么在《R for Dummies》的第60页上有解释。您还可以查看?na.omit，或者查看SO上的这个问题：https://dev59.com/cG445IYBdhLWcg3wcJ2O - Joris Meys
显示剩余4条评论
43

只需将以下内容添加到您的程序中，即可获得曲线下的面积：

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

来自?trapz:

使用梯形法规则和基准点x对函数进行积分的近似方法完全匹配。

2
细节总是受欢迎的，特别是当答案已经被接受时。 - Nikana Reklawyks
1
请注意，如果您的 x 值是递减的，则 trapz() 将给出负值。请参考 x<-1:10 与 x<-10:1 的区别。 - Matt
24

还有三个选项，其中一个使用样条方法，另一个使用辛普森规则...

# 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/one-dimensional-integrals/

1
有另一个叫做“flux”的包。它具有与“MESS”相同的函数名称，“auc()”。值得一试！ - Fábio
23

好的，所以我有点晚到派对现场，但是在查看答案时发现缺少一个简单的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)。

1
我总是更喜欢简单的 R 解决方案 :) - Verbal
6
在药代动力学（PK）领域，计算不同类型的AUC是一项常见而基础的任务。药代动力学中有许多不同的AUC计算方法，例如：
  • AUC0-t = 从零到时间t的累积面积
  • AUC0-last = 从零到最后一个时间点的累积面积（可能与上面相同）
  • AUC0-inf = 从零到无限时间的累积面积
  • AUCint = 在时间间隔内的累积面积
  • AUCall = 数据存在的整个时间段内的累积面积
其中一个最好的可以进行这些计算的软件包是较新的Pfizer公司的PKNCA软件包。快去试试吧！
1

Joris Meys的回答很好，但我在从我的样本中删除NAs方面遇到了困难。这里是我编写的小函数来处理它们：

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)) 希望这能帮助像我这样的新手 :-)
编辑：添加了边界
-6
您可以使用ROCR软件包，以下代码将为您提供AUC值：
pred <- prediction(classifier.labels, actual.labs)
attributes(performance(pred, 'auc'))$y.values[[1]]
7
OP不想计算ROC曲线和其AUC，而是任意曲线下的面积。 - Calimo
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接