gpt4 book ai didi

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

转载 作者:行者123 更新时间:2023-12-02 01:24:54 30 4
gpt4 key购买 nike

我想计算由一组实验值定义的曲线下面积。我创建了一个函数来使用辛普森规则计算 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)

最佳答案

处理偶数个点并仍然实现精度的一种推荐方法是将辛普森 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

关于r - 使用 R 中的辛普森规则计算曲线下面积,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/74977397/

30 4 0
Copyright 2021 - 2024 cfsdn All Rights Reserved 蜀ICP备2022000587号
广告合作:1813099741@qq.com 6ren.com