67

我想计算曲线下的面积以进行积分而不定义诸如 in 之类的函数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)查看波动性微笑。有没有办法整合这个绘制的“曲线”？

4

7 回答 7

49

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

于 2011-02-10T09:07:10.940 回答
37

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

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

来自?trapz：

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

于 2012-03-01T21:08:38.900 回答
22

另外三个选项，包括一个使用样条方法和一个使用辛普森规则......

# 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/

于 2013-01-29T21:34:22.410 回答
21

好的，所以我在聚会上来的有点晚，但是仔细检查答案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)如果这更有意义，您可以使用。

于 2015-05-16T21:19:22.250 回答
5

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

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

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

于 2016-05-03T20:42:07.563 回答
1

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

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

编辑：添加界限

于 2017-06-15T09:53:02.603 回答
-5

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

pred <- prediction(classifier.labels, actual.labs)
attributes(performance(pred, 'auc'))$y.values[[1]]
于 2012-10-26T08:26:16.067 回答