r - 使用 R 中的辛普森规则计算曲线下面积

Question

我想计算由一组实验值定义的曲线下面积。我创建了一个函数来使用辛普森规则计算 AUC 的近似值，如我在 this post 中看到的那样。但是，该函数仅在接收奇数长度的向量时才起作用。当输入向量的长度为偶数时，如何修改代码以添加最后一个梯形的面积。

AUC <- function(x, h=1){
  # AUC function computes the Area Under the Curve of a time serie using
  # the Simpson's Rule (numerical method).
  # https://link.springer.com/chapter/10.1007/978-1-4612-4974-0_26
  
  # Arguments
  # x: (vector) time serie values
  # h: (int) temporal resolution of the time serie. default h=1
  
  n = length(x)-1
  
  xValues = seq(from=1, to=n, by=2)
  
  sum <- list()
  for(i in 1:length(xValues)){
    n_sub <- xValues[[i]]-1
    n <- xValues[[i]]
    n_add <- xValues[[i]]+1
    
    v1 <- x[[n_sub+1]]
    v2 <- x[[n+1]]
    v3 <- x[[n_add+1]]
    
    s <- (h/3)*(v1+4*v2+v3)
    sum <- append(sum, s)
  }
  sum <- unlist(sum)
  
  auc <- sum(sum)
  
  return(auc)
  
}

这里是一个数据示例:

smoothed = c(0.3,0.317,0.379,0.452,0.519,0.573,0.61,0.629,0.628,0.613,0.587,0.556,0.521,
             0.485,0.448,0.411,0.363,0.317,0.273,0.227,0.185,0.148,0.12,0.103,0.093,0.086,
             0.082,0.079,0.076,0.071,0.066,0.059,0.053,0.051,0.052,0.057,0.067,0.081,0.103,
             0.129,0.165,0.209,0.252,0.292,0.328,0.363,0.398,0.431,0.459,0.479,0.491,0.494,
             0.488,0.475,0.457,0.43,0.397,0.357,0.316,0.285,0.254,0.227,0.206,0.189,0.181,
             0.171,0.157,0.151,0.162,0.192,0.239)

Answer 1

处理偶数个点并仍然实现精度的一种推荐方法是将辛普森 1/3 规则与辛普森 3/8 规则结合起来，后者可以处理偶数个点。这些方法可以在(至少一本或多本)关于数值方法的工程教科书中找到。

但是，实际上，您可以编写一个代码块来检查数据长度并在末尾添加一个梯形，正如您链接到的帖子的最后一条评论中所建议的那样。我不认为它一定像结合辛普森的 1/3 和 3/8 规则一样精确，但对于许多应用程序来说它可能是合理的。

我会仔细检查下面的代码编辑，但这是基本想法。

AUC <- function(x, h=1){
  # AUC function computes the Area Under the Curve of a time serie using
  # the Simpson's Rule (numerical method).
  # https://link.springer.com/chapter/10.1007/978-1-4612-4974-0_26
  
  # Arguments
  # x: (vector) time serie values
  # h: (int) temporal resolution of the time serie. default h=1
  
  #jh edit: check for even data length
  #and chop off last data point if even
  nn = length(x)
  if(length(x) %% 2 == 0){
  xlast = x[length(x)]
  x = x[-length(x)]
  }

  n = length(x)-1
  
  xValues = seq(from=1, to=n, by=2)
  
  sum <- list()
  for(i in 1:length(xValues)){
    n_sub <- xValues[[i]]-1
    n <- xValues[[i]]
    n_add <- xValues[[i]]+1
    
    v1 <- x[[n_sub+1]]
    v2 <- x[[n+1]]
    v3 <- x[[n_add+1]]
    
    s <- (h/3)*(v1+4*v2+v3)
    sum <- append(sum, s)
  }
  sum <- unlist(sum)
  
  auc <- sum(sum)
  ##jh edit: add trapezoid for last two data points to result
  if(nn %% 2 == 0){
  auc <- auc + (x[length(x)] + xlast)/2 * h
  }
  
  
  return(auc)
  
}


sm = smoothed[-length(smoothed)]
length(sm)
[1] 70
#even data as an example
AUC(sm)
[1] 20.17633
#original odd data
AUC(smoothed)
[1] 20.389