计算曲线下的面积

我想计算曲线下的面积来完成整合，而不需要像在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

通过观察大量的梯形图，每次在x_i ， x_{i+1} ， y{i+1}和y_i之间进行界限，AUC可以非常容易地估计出来。 使用动物园包的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() ）

关于你的后续工作：如果你没有正式的function，你会如何绘制它？ 所以如果你只有价值观，你唯一可以接近的就是一个确定的积分。 即使在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 （使用辛普森规则）应该是首选。 这些DIY版本（以及使用正交规则的附加方法）在这里： http : //www.r-bloggers.com/one-dimensional-integrals/

Answer 4

好的，所以我在聚会上迟到了一些，但是回答一个简单的问题解决scheme就没有了。 在这里，简单而干净：

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

OP的解决scheme读取为：

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

这通过采用“左”和“右”y值的平均值，使用梯形法有效地计算面积。

注意：正如@Joris已经指出的那样，如果这样做更有意义，你可以使用abs(y) 。

Answer 5

在药代动力学（PK）世界中，计算不同types的AUC是一个常见的基本任务。 药代动力学有很多不同的AUC计算方法，如

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

做这些计算的最好的软件包之一是来自辉瑞公司的人员的相对较新的软件包PKNCA 。 一探究竟。

Answer 6

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))上： myDF$auc = apply(myDF, MARGIN=1, FUN=auc, x=c(0,5,10,15,20))

希望它可以帮助像我这样的noobs 🙂

编辑：增加了界限

Answer 7

您可以使用ROCR软件包，其中以下行将给你AUC：

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