如何计算正弦曲线两端的面积

给定这个数据集：

y<-c(-13,16,35,40,28,36,43,33,40,33,22,-5,-27,-31,-29,-25,-26,-31,-26,-24,-25,-29,-23,4)
t<-1:24

我的目标是计算两个面积。第一个区域将仅集成来自零线上方的曲线第一部分的数据。第二个区域将整合来自零线下方的曲线第二部分的数据。

首先，我想对该数据拟合正弦波。使用这个优秀的答案：

https://stats.stackexchange.com/questions/60994/fit-a-sinusoidal-term-to-data

我能够拟合正弦波（我将使用带有二次谐波的周期波，它看起来更适合）

ssp <- spectrum(y)  
per <- 1/ssp$freq[ssp$spec==max(ssp$spec)]
reslm <- lm(y ~ sin(2*pi/per*t)+cos(2*pi/per*t))
summary(reslm)

rg <- diff(range(y))
plot(y~t,ylim=c(min(y)-0.1*rg,max(y)+0.1*rg))
lines(fitted(reslm)~t,col=4,lty=2)   # dashed blue line is sin fit

# including 2nd harmonic really improves the fit
reslm2 <- lm(y ~ sin(2*pi/per*t)+cos(2*pi/per*t)+sin(4*pi/per*t)+cos(4*pi/per*t))
summary(reslm2)
lines(fitted(reslm2)~t,col=3)    # solid green line is periodic with second harmonic
abline(h=0,lty=2)

接下来，我想计算仅为正的曲线下的面积，以及仅为负的曲线下的面积。我很幸运地使用 Bolstad2 和 Mess 包中的 AUC 函数查看了类似的答案。但我的数据点并没有整齐地落在零线上，而且我不知道如何将正弦函数分解为仅在零线上方和仅在零线下方的区域。

首先要事。为了获得精确的计算，您需要使用二次谐波傅里叶的精确函数。其次，谐波函数的美妙之处在于它们是重复的。因此，如果你想找到你的函数在哪里达到 0，你只需要把你的区间扩大到，这样你就一定能找到 2 个以上的根。

首先我们从回归模型中得到精确的函数

fourierfnct <- function(t){
  fnct <- reslm2$coeff[1]+
    reslm2$coeff[2]*sin(2*pi/per*t)+
    reslm2$coeff[3]*cos(2*pi/per*t)+
    reslm2$coeff[4]*sin(4*pi/per*t)+
    reslm2$coeff[5]*cos(4*pi/per*t)
  return(fnct)
}

其次，你可以编写一个可以求根的函数（其中函数为0）。 R 提供了一个 uniroot 函数，您可以使用它来查找循环中的多个根。

manyroots <- function(f,inter,period){
  roots <- array(NA, inter)
  for(i in 1:(length(inter)-1)){
    roots[i] <- tryCatch({
      return_value <- uniroot(f,c(inter[i],inter[i+1]))$root
    }, error = function(err) {
      return_value <- -1
    })
  }
  retroots <- roots[-which(roots==-1)]
  return(retroots)
}

然后，您只需计算根，并使用它们来跨这些边界对函数进行积分。

roots <- manyroots(fourierfnct,seq(0,25),per)
integrate(fourierfnct, roots[1],roots[2])
#300.6378 with absolute error < 3.3e-12
integrate(fourierfnct, roots[2],roots[3])
#-284.6378 with absolute error < 3.2e-12