计算曲线下的面积

Posted 2023-03-12

【中文标题】计算曲线下的面积

【英文标题】：Calculate the Area under a Curve

【发布时间】：2011-06-24 16:19:54

【问题描述】：

我想以进行积分而不定义函数，例如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) 来查看波动微笑。有没有办法整合这条绘制的“曲线”？

【问题讨论】：

看看这个相关问题：***.com/questions/4903092/calculate-auc-in-r

@Andrie : 这是不同类型的 AUC...

【参考方案1】：

AUC 很容易通过查看大量梯形图来近似，每次都在x_i、x_i+1、yi+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() 计算定积分。仅当您也可以定义形式函数时，才能绘制形式函数。

【讨论】：

谢谢，这行得通。还有一种方法可以绘制积分吗？我的意思是，如果我有一条像波动率微笑这样的曲线，我应该能够绘制它的积分，这也是一条曲线。

这是测试我的 pdf 总和为 1 的好方法。谢谢！

这很好，但如果缺少某些值，公式将不再有效。

@DanChaltiel 如果缺少某些值，则无法知道曲线下的实际面积是多少。所以这对我来说似乎不是问题。如果您想忽略缺失的数据，只需在计算前删除缺失的观察值即可。

@JorisMeys 如果您有 10 个 x 值而只有 9 个 y 值，那么如果您不计算缺失值，您可以得到一个非常好的 AUC 近似值。删除所有只有一个 NA 的样本对我来说似乎是一种浪费。

【参考方案2】：

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

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

来自?trapz：

这种方法与积分的近似值完全匹配 使用带有基点 x 的梯形规则的函数。

【讨论】：

欢迎提供详细信息，尤其是当答案已被接受时。

请注意，如果您的 x 值正在下降，trapz() 会给您一个负值。请参阅x<-1:10 与 x<-10:1。

【参考方案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/

【讨论】：

还有一个名为“flux”的包。它具有与“MESS”相同的函数名称“auc()”。值得一试！

【参考方案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)。

【讨论】：

我总是更喜欢简单的R 解决方案:)

【参考方案5】：

在药代动力学 (PK) 领域，计算不同类型的 AUC 是一项常见且基本的任务。药代动力学有很多不同的 AUC 计算，例如

AUC0-t = AUC 从零到时间 t AUC0-last = AUC 从零到最后一个时间点（可能同上） AUC0-inf = AUC 从零到时间无穷 AUCint = 一段时间内的 AUC AUCall = 存在数据的整个时间段内的 AUC

进行这些计算的最佳软件包之一是辉瑞公司提供的相对较新的软件包PKNCA。看看吧。

【参考方案6】：

Joris Meys's answer 很棒，但我很难从样本中删除 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))

希望它可以帮助像我这样的菜鸟:-)

编辑：添加边界

【参考方案7】：

您可以使用 ROCR 包，其中以下几行将为您提供 AUC：

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

【讨论】：

OP不想计算ROC曲线及其AUC，而是计算任意曲线下的面积。