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

Question

鉴于此数据集：

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 函数查看了类似的答案。但是我的数据点并没有整齐地落在零线上，我不知道如何将正弦函数分解成仅在零线上方和零线下方的区域。

Answer 1

这可能不是您正在寻找的解决方案，但您可以试试这个：

# Create a new t vector but with more subdivisions
t2 = seq(1,24,length.out = 10000)

# Evaluate your model on this t2
y2 = predict(reslm2, newdata = data.frame(t = t2))

lines(t2[y2>=0],y2[y2>=0],col="red")

# Estimate the area where the curve is greater than 0 
sum(diff(t2)[1]*y2[y2>0])

# Estimate the area where the curve is less than 0
sum(diff(t2)[1]*y2[y2<0])

Answer 2

要事第一。要获得精确的计算，您需要使用二次谐波傅里叶的精确函数。其次，谐波函数的优点在于它们是重复的。所以如果你想找到你的函数到达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

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

How to calculate the area under each end of a sine curve

r

area

trigonometry