在 R 中使用 lm、nls 和 nls2 进行正弦曲线拟合答案

【问题标题】：Sine curve fit using lm, nls and nls2 in R在 R 中使用 lm、nls 和 nls2 进行正弦曲线拟合
【发布时间】：2021-09-24 17:56:35
【问题描述】：

我获得了数据（运动适应 =y 与延迟 =t 的关系），我希望它看起来像一个正弦波。我正在尝试：

在我的数据中拟合正弦曲线
估计最佳模型/参数。

我已经阅读了几篇 here、here 和 here 的帖子，但我仍在苦苦挣扎。

1) 使用 lm

代码：

t<-c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9)
y<-c(0.310 ,0.630 ,0.430 ,0.245, 0.650 ,0.085 ,0.370, 0.560 ,0.250, 0.520)
reslm <- lm(y ~ sin(pi/2*t)+ cos(pi/2*t)) #my period is supposed to be 4, so period equals to pi/2
summary(reslm)
rg<-(max(y)-min(y)/2)
plot(y~t)
lines(fitted(reslm)~t,col=4,lty=2)

输出：

lm(formula = y ~ sin(pi/2 * t) + cos(pi/2 * t))

Residuals:
     Min       1Q   Median       3Q      Max 
-0.32450 -0.13956 -0.00325  0.14819  0.24450 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)   0.404375   0.067993   5.947 0.000572 ***
sin(pi/2 * t) 0.005125   0.095190   0.054 0.958567    
cos(pi/2 * t) 0.001125   0.095190   0.012 0.990900    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2107 on 7 degrees of freedom
Multiple R-squared:  0.0004303, Adjusted R-squared:  -0.2852 
F-statistic: 0.001507 on 2 and 7 DF,  p-value: 0.998

图表：

问题：

我很困惑，如何改变幅度和相移？
如何使用这种方法提高我的合身度？

2) 使用 nls

我使用了公式y(t) = A*sin(Omega*t + Phi) + C，其中A 是幅度，Omega 是周期，Phi 是相移，C 是中线。

代码：

t<-c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9)
y<-c(0.310 ,0.630 ,0.430 ,0.245, 0.650 ,0.085 ,0.370, 0.560 ,0.250, 0.520)

A<- (max(y)-min(y)/2)
C<-((max(y)+min(y))/2)

res1<- nls(y ~ A*sin(omega*t+phi)+C, data=data.frame(t,y), start=list(A=A,omega=pi/2,phi=0,C=C))
summary(res1)
co <- coef(res1)
resid(res1)
sum(resid(res1)^2)
fit <- function(x, a, b, c, d) {a*sin(b*x+c)+d}
# Plot result
plot(x=t, y=y)
curve(fit(x, a=co["A"], b=co["omega"], c=co["phi"], d=co["C"]), add=TRUE ,lwd=2, col="steelblue")

输出：

Formula: y ~ A * sin(omega * t + phi) + C

Parameters:
       Estimate Std. Error t value Pr(>|t|)    
A       0.21956    0.03982   5.513   0.0015 ** 
omega   2.28525    0.07410  30.841 7.72e-08 ***
phi   -32.57364    0.40375 -80.678 2.44e-10 ***
C       0.41146    0.02926  14.061 8.07e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.09145 on 6 degrees of freedom

Number of iterations to convergence: 18 
Achieved convergence tolerance: 9.705e-06

图表：

问题：

这种方法似乎效果更好一些，但是如何使用这种方法提高拟合度？我尝试手动更改一些参数，例如，更改相移 (phi) 没有太大作用或导致错误（请参阅第 3 部分）。

3) 使用 nls 和 nsl2 来调整我的模型

代码：

###nls2
t<-c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9)
y<-c(0.310 ,0.630 ,0.430 ,0.245, 0.650 ,0.085 ,0.370, 0.560 ,0.250, 0.520)
A <- (max(y)-min(y)/2)
C<-((max(y)+min(y))/2)
pp <- expand.grid(omega=(c(2.094395, 1.570796,  1.256637)), phi=(-1:1), A=A, C=C) # omega = 2*pi/3, pi/2 , 2*pi/5
View(pp)
pp1<-data.frame(pp)
res2<- nls2(y ~ A*sin(omega*t+phi)+C, data=data.frame(t,y), start=pp1,  algorithm = "brute-force")
res2
summary(res2)
co <- coef(res2)
resid(res2)
sum(resid(res2)^2)
fit <- function(x, a, b, c, d) {a*sin(b*x+c)+d}
# Plot result 
plot(x=t, y=y)
curve(fit(x, a=co["A"], b=co["omega"], c=co["phi"], d=co["C"]), add=TRUE ,lwd=2, col="steelblue")

#optimisation
res3<-nls2(y ~ A*sin(omega*t+phi)+C, start = res2)
res3
summary(res3)
co3 <- coef(res3)
resid(res3)
sum(resid(res3)^2)
fit <- function(x, a, b, c, d) {a*sin(b*x+c)+d}
# Plot result
plot(x=t, y=y)
curve(fit(x, a=co3["A"], b=co["omega"], c=co3["phi"], d=co3["C"]), add=TRUE ,lwd=2, col="steelblue")

输出：

第一次尝试（nls2 model1）：

model: y ~ A * sin(omega * t + phi) + C
   data: data.frame(t, y)
 omega    phi      A      C 
2.0944 0.0000 0.6075 0.3675 
 residual sum-of-squares: 0.8545

Number of iterations to convergence: 9 
Achieved convergence tolerance: NA
> summary(res2)

Formula: y ~ A * sin(omega * t + phi) + C

Parameters:
      Estimate Std. Error t value Pr(>|t|)    
omega  2.09440    0.08453  24.776 2.84e-07 ***
phi    0.00000    0.46494   0.000   1.0000    
A      0.60750    0.17851   3.403   0.0144 *  
C      0.36750    0.12044   3.051   0.0225 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3774 on 6 degrees of freedom

Number of iterations to convergence: 9 
Achieved convergence tolerance: NA

图表：

第二次尝试（nls2 model2）：

Nonlinear regression model
  model: y ~ A * sin(omega * t + phi) + C
   data: <environment>
  omega     phi       A       C 
 2.2852 -1.1577  0.2196  0.4115 
 residual sum-of-squares: 0.05018

Number of iterations to convergence: 12 
Achieved convergence tolerance: 8.075e-06
> summary(res3)

Formula: y ~ A * sin(omega * t + phi) + C

Parameters:
      Estimate Std. Error t value Pr(>|t|)    
omega  2.28524    0.07410  30.841 7.72e-08 ***
phi   -1.15769    0.40375  -2.867   0.0285 *  
A      0.21956    0.03982   5.513   0.0015 ** 
C      0.41146    0.02926  14.061 8.07e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.09145 on 6 degrees of freedom

Number of iterations to convergence: 12 
Achieved convergence tolerance: 8.075e-06

图表：

问题：

因此，我似乎误解了 nls2 所做的事情，因为我发现与第 2 部分完全相同的结果。我仍然不知道哪些参数是最好的，我该怎么做？

4) 使用 nls，通过循环几个参数来调整我的模型

代码：

t<-c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9)
y<-c(0.310 ,0.630 ,0.430 ,0.245, 0.650 ,0.085 ,0.370, 0.560 ,0.250, 0.520)
A <- (max(y)-min(y)/2)
C<-((max(y)+min(y))/2)
pp <- expand.grid(omega=(c(2.094395, 1.570796,  1.256637)), phi=(-1:1), A=A, C=C) # omega = 2*pi/3, pi/2 , 2*pi/5
#View(pp)

fit_AIC<- vector()
fit_BIC<- vector()

coef_A<- vector()
coef_ome<- vector()
coef_phi<- vector()
coef_C<- vector()
RSS<-vector()

for (ii in 1:nrow(pp))
{
    res<- nls(y ~ A*sin(omega*t+phi)+C, data=data.frame(t,y), start=list(A=pp$A[ii],omega=pp$omega[ii],phi=pp$phi[ii],C=pp$C[ii]), trace = TRUE)
    fit_AIC[ii]<-AIC(res)
    fit_BIC[ii]<-BIC(res)

    coef_A[ii]<- coef(res)[1]
    coef_ome[ii]<- coef(res)[2]
    coef_phi[ii]<- coef(res)[3]
    coef_C[ii]<- coef(res)[4]
    RSS<-sum(resid(res)^2)

}

results<-data.frame(RSS, fit_AIC,  fit_BIC,  coef_A,  coef_ome,  coef_phi,  coef_C)
View(results)

输出：

我收到此错误：

1.405742 :   0.607500  2.094395 -1.000000  0.367500
0.1448148 :   0.1563179  2.1441802 -0.9937729  0.4172079

...


0.05018035 :   0.2195573  2.2852482 -1.1577097  0.4114573
2.085664 :  0.607500 1.570796 1.000000 0.367500
0.3104012 :  0.01321257 1.60518024 0.83201816 0.40437498
0.3098916 :   0.0180852  3.0888764 -5.9933691  0.4060743
Error in nls(y ~ A * sin(omega * t + phi) + C, data = data.frame(t, y),  : 
  le pas 0.000488281 est devenu inférieur à 'minFactor' de 0.000976562


        RSS    fit_AIC    fit_BIC     coef_A  coef_ome    coef_phi    coef_C
1 0.05018035 -14.568398 -13.055473 0.21955754 2.2852455  -1.1576955 0.4114573
2 0.05018035   2.753153   4.266079 0.07487110 0.8575642   0.2299909 0.3916769
3 0.05018035   2.753153   4.266079 0.07487109 0.8575736   0.2299951 0.3916763
4 0.05018035 -14.568398 -13.055473 0.21955763 2.2852443  -1.1576894 0.4114573
5 0.05018035 -14.568398 -13.055473 0.21955729 2.2852490 -32.5736406 0.4114573
6 0.05018035   2.753153   4.266079 0.07487105 0.8575619   0.2300021 0.3916770
7 0.05018035 -14.568398 -13.055473 0.21955735 2.2852482  -1.1577097 0.4114573

问题：

所以这个错误似乎是因为我的初始参数错误。这个对吗？但是，如果大多数参数都不起作用，我该如何估计最佳参数？
另外，我不明白为什么尽管参数不同，RSS 总是相同的
为什么我只观察到 2 个不同的 AIC 和 2 个不同的 BIC 而模型不同？

任何形式的帮助将不胜感激，谢谢。

【问题讨论】：

您在最后表的所有行中最小化相同的函数，并且该函数具有多个局部最小值。例如，看这个图：curve(x^4 - 4 * x^3 - 2 * x^2 + 9 * x + 9, -2, 5)。如果您使用靠近图表左侧的起始值最小化该函数，您可能会得到第一个局部最小值，如果您从靠近右侧的起始值开始，您可能会得到第二个。
您好，感谢您的回复，但我有点困惑。你在这里回答我的哪个问题？我应该怎么做才能找到最佳参数？还是我应该接受第 2 点的结果作为我的最佳选择？谢谢

标签： r trigonometry lm nls

【解决方案1】：

在第一个模型上，您正在有效地拟合模型 y = A*sin(pi/2 * t) + B*cos(pi/2 * t) + C，并且您无法使用该方法更改相移。在接下来的步骤中，您可以拟合一个适当的方程。 nls 使用提供的模型执行非线性回归，并且 nls2 通过创建初始参数网格并执行多个 nls 调用添加了一个额外的步骤。在最后一种方法中，您编写了与 nls2 类似的策略。但是一些初始条件会导致错误。

您最小化的目标函数具有多个局部最小值，就像示例中的函数具有多个具有不同坐标的最小值一样。

此外，当您将等式更改为 y ~ A * sin(omega * t + phi) + C 时，您需要确定四个参数，并且只有九个数据点需要添加更多数据点才能更好地估计参数。您还可以将数据标准化为 0-1 范围。

希望对你有帮助

【讨论】：

【解决方案2】：

使用 lm：您的经期约为 2.75 而不是 4。

  t<-c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9)
  y<-c(0.310 ,0.630 ,0.430 ,0.245, 0.650 ,0.085 ,0.370, 0.560 ,0.250, 0.520)

  reslm <- lm(y ~ cos(2* pi/2.75*t) + sin(2* pi/2.75*t))

  summary(reslm)

  coef1 <- coef(reslm)

  plot(y~t)

  t1 <- seq(0,9,.1)
  y1 <- coef1[2] * cos(2*pi/2.75*t1) + coef1[3] * sin(2* pi/2.75*t1) + coef1[1]

  lines(t1,y1,col=2)

【讨论】：