机器算法验证 - 用 R 中的 lm 拟合正弦波以获取昼夜节律活动频率？ - 吾爱随笔录

用 R 中的 lm 拟合正弦波以获取昼夜节律活动频率？

机器算法验证 r 回归时间序列循环统计余弦

2022-03-04 01:37:27

我正在尝试在一些活动数据上拟合正弦波，就像这篇文章一样。我已经设法为一种情况获得了一个合理的图表：

然而，当我绘制其他条件时，波浪看起来非常奇怪，而不是它应该看起来的样子——它应该像第一张图像一样漂亮和平滑。此外，它应该像第一张图像一样有两个峰和一个谷。

这里出了什么问题？我能做些什么来解决这个问题吗？第二个条件具有相同的 48 小时，依此类推。

我的数据在 48 小时内（连续两天，因此每天应该有一个高峰和一个低谷）具有每小时测量值。代码从这里改编。

b$zt <- b$zt/48
lmfit <- lm(data = b,
            walking ~ sin(2*pi*b$zt) + cos(2*pi*b$zt))
b0 <- coef(lmfit)[1]
alpha <- coef(lmfit)[2]
beta <- coef(lmfit)[3]

pframe <- data.frame(zt=seq(min(b$zt),max(b$zt),length=612))
pframe$walking <- predict(lmfit,newdata=pframe)
library(ggplot2); theme_set(theme_bw())
ggplot(b,aes(zt,walking))+
  geom_point(alpha=0.2) + geom_line(data=pframe,colour="red")

以下是数据集中的相关部分：

 dput(new$zt)
c(6L, 6L, 6L, 6L, 7L, 7L, 7L, 7L, 8L, 8L, 8L, 8L, 9L, 9L, 9L, 
9L, 11L, 11L, 11L, 11L, 12L, 12L, 12L, 12L, 15L, 15L, 15L, 15L, 
17L, 17L, 17L, 17L, 23L, 23L, 23L, 23L, 24L, 24L, 24L, 24L, 1L, 
1L, 1L, 1L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 7L, 7L, 7L, 7L, 8L, 
8L, 8L, 8L, 9L, 9L, 9L, 9L, 10L, 10L, 10L, 10L, 11L, 11L, 11L, 
11L, 12L, 12L, 12L, 12L, 13L, 13L, 13L, 13L, 14L, 14L, 14L, 14L, 
15L, 15L, 15L, 15L, 17L, 17L, 17L, 17L, 19L, 19L, 19L, 19L, 24L, 
24L, 24L, 24L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 7L, 7L, 7L, 7L, 
8L, 8L, 8L, 8L, 9L, 9L, 9L, 9L, 10L, 10L, 10L, 10L, 11L, 11L, 
11L, 11L, 12L, 12L, 12L, 12L, 22L, 22L, 22L, 22L, 24L, 24L, 24L, 
24L, 2L, 2L, 2L, 2L, 5L, 5L, 5L, 6L, 6L, 6L, 7L, 7L, 7L, 8L, 
8L, 8L, 9L, 9L, 9L, 10L, 10L, 10L, 11L, 11L, 11L, 12L, 12L, 12L, 
14L, 14L, 14L, 15L, 15L, 15L, 16L, 16L, 16L, 17L, 17L, 17L, 18L, 
18L, 18L, 19L, 19L, 19L, 21L, 21L, 21L, 22L, 22L, 22L, 23L, 23L, 
23L, 24L, 24L, 24L, 2L, 2L, 2L, 31L, 31L, 31L, 31L, 32L, 32L, 
32L, 32L, 33L, 33L, 33L, 33L, 35L, 35L, 35L, 35L, 36L, 36L, 36L, 
36L, 39L, 39L, 39L, 39L, 41L, 41L, 41L, 41L, 42L, 42L, 42L, 42L, 
44L, 44L, 44L, 44L, 45L, 45L, 45L, 45L, 46L, 46L, 46L, 46L, 25L, 
25L, 25L, 25L, 26L, 26L, 26L, 26L, 28L, 28L, 28L, 28L, 30L, 30L, 
30L, 30L, 35L, 35L, 35L, 35L, 29L, 29L, 29L, 29L, 30L, 30L, 30L, 
30L, 31L, 31L, 31L, 31L, 32L, 32L, 32L, 32L, 33L, 33L, 33L, 33L, 
34L, 34L, 34L, 34L, 35L, 35L, 35L, 35L, 36L, 36L, 36L, 36L, 38L, 
38L, 38L, 38L, 41L, 41L, 41L, 41L, 42L, 42L, 42L, 42L, 45L, 45L, 
45L, 45L, 46L, 46L, 46L, 46L, 29L, 29L, 29L, 30L, 30L, 30L, 31L, 
31L, 31L, 32L, 32L, 32L, 33L, 33L, 33L, 34L, 34L, 34L, 35L, 35L, 
35L, 36L, 36L, 36L, 37L, 37L, 37L, 38L, 38L, 38L, 39L, 39L, 39L, 
40L, 40L, 40L, 41L, 41L, 41L, 42L, 42L, 42L, 44L, 44L, 44L, 45L, 
45L, 45L, 25L, 25L, 25L, 13L, 13L, 13L, 13L, 14L, 14L, 14L, 14L, 
16L, 16L, 16L, 16L, 18L, 18L, 18L, 18L, 19L, 19L, 19L, 19L, 20L, 
20L, 20L, 20L, 21L, 21L, 21L, 21L, 22L, 22L, 22L, 22L, 13L, 13L, 
13L, 13L, 14L, 14L, 14L, 14L, 15L, 15L, 15L, 15L, 17L, 17L, 17L, 
17L, 18L, 18L, 18L, 18L, 21L, 21L, 21L, 21L, 16L, 16L, 16L, 16L, 
18L, 18L, 18L, 18L, 20L, 20L, 20L, 20L, 21L, 21L, 21L, 21L, 22L, 
22L, 22L, 22L, 23L, 23L, 23L, 23L, 13L, 13L, 13L, 13L, 14L, 14L, 
14L, 14L, 15L, 15L, 15L, 15L, 16L, 16L, 16L, 16L, 17L, 17L, 17L, 
17L, 18L, 18L, 18L, 18L, 19L, 19L, 19L, 19L, 20L, 20L, 20L, 20L, 
21L, 21L, 21L, 21L, 14L, 14L, 14L, 14L, 16L, 16L, 16L, 16L, 17L, 
17L, 17L, 17L, 18L, 18L, 18L, 18L, 21L, 21L, 21L, 21L, 22L, 22L, 
22L, 22L, 13L, 13L, 13L, 20L, 20L, 20L, 13L, 13L, 13L, 14L, 14L, 
14L, 17L, 17L, 17L, 21L, 21L, 21L)

> dput(new$walking)
c(8, 1.166666667, 4.666666667, 2, 4.833333333, 2, 1.833333333, 
4, 2.5, 1.333333333, 4.166666667, 1.833333333, 7, 1.166666667, 
9.5, 1.166666667, 6.833333333, 5.166666667, 4.333333333, 2.333333333, 
5.166666667, 2.333333333, 2, 1.833333333, 0, 0, 0, 0, 0, 0, 1.5, 
0, 0, 0, 1, 0.333333333, 0.333333333, 1.5, 2.5, 2.333333333, 
4.333333333, 8.833333333, 6.666666667, 6.166666667, 2.166666667, 
0.333333333, 0, 0, 9.666666667, 12.5, 13.5, 16.16666667, 0, 2.166666667, 
0.666666667, 0.5, 0, 2.333333333, 0, 0, 0, 1.5, 3, 0.833333333, 
3, 7.5, 5, 3.166666667, 7.666666667, 7.833333333, 2.666666667, 
2.5, 2.333333333, 1, 1.666666667, 0.333333333, 1.5, 0.666666667, 
0.333333333, 3.166666667, 1.833333333, 0, 0, 0, 0, 1, 0, 0, 0, 
1.166666667, 0, 0, 0, 0.833333333, 0.666666667, 0.666666667, 
0, 4, 4.5, 3.833333333, 0, 3, 0.333333333, 0, 4.5, 1.333333333, 
5, 3, 0, 0.666666667, 0.333333333, 1.5, 3, 5.5, 1, 6.166666667, 
1.5, 3.333333333, 15.33333333, 4.166666667, 2, 1.333333333, 7.166666667, 
1.5, 1.333333333, 5.333333333, 4.666666667, 6.166666667, 0.333333333, 
1.166666667, 0.833333333, 1.5, 0, 0, 0, 0.166666667, 1.166666667, 
2, 1.333333333, 4.333333333, 3.833333333, 1.333333333, 0.333333333, 
0.666666667, 1.5, 4.833333333, 0.666666667, 0.666666667, 0, 4, 
1.5, 0.833333333, 1.166666667, 2.166666667, 0.5, 0.666666667, 
0, 0, 0.833333333, 2.5, 1, 0, 0.833333333, 3.5, 2.166666667, 
1, 1.333333333, 0.5, 0.333333333, 0.166666667, 0, 0, 0.166666667, 
0.333333333, 0.333333333, 0.166666667, 0.5, 0, 0.166666667, 0, 
0, 0, 0, 0, 0.166666667, 0, 0, 0, 0.166666667, 0, 0.166666667, 
0, 0, 0, 0.5, 1.333333333, 4, 2.666666667, 3.833333333, 2.5, 
10, 1.166666667, 5, 2.5, 2.333333333, 3.833333333, 4.833333333, 
2, 2.833333333, 0, 0.666666667, 0, 0, 3, 3.5, 5.166666667, 3.166666667, 
3.5, 2.833333333, 2.166666667, 8.166666667, 0.166666667, 0.5, 
0, 0, 0.333333333, 0.166666667, 0.333333333, 0, 0, 0.666666667, 
0, 0, 0.166666667, 0.333333333, 0.666666667, 0, 0.166666667, 
0.166666667, 1.166666667, 0.666666667, 3, 6.5, 4, 6.166666667, 
0, 3.166666667, 1.333333333, 2, 0, 0, 1.666666667, 0, 1.833333333, 
0, 0.333333333, 0, 0, 0.166666667, 0.5, 0, 3.666666667, 1.166666667, 
3.833333333, 2.5, 0, 0, 0, 0, 1.666666667, 1.166666667, 0.666666667, 
0, 0.833333333, 2.833333333, 1.333333333, 0.5, 2.166666667, 1.5, 
1.666666667, 4.333333333, 12.16666667, 9.833333333, 3.666666667, 
3.333333333, 1.166666667, 6.666666667, 1.833333333, 1.166666667, 
5.5, 4.166666667, 5.833333333, 16.83333333, 2.602230483, 6.319702602, 
1.672862454, 4.089219331, 0, 0.666666667, 0, 0, 0.166666667, 
0.166666667, 0, 0, 0, 0, 0, 1, 0, 0.166666667, 0.333333333, 0.833333333, 
18.33333333, 5.833333333, 1.5, 4.5, 3.5, 7.166666667, 2.5, 7.833333333, 
3.666666667, 6.666666667, 0.333333333, 5.666666667, 5.5, 2.166666667, 
5.333333333, 2.333333333, 1.166666667, 1.166666667, 0, 0.833333333, 
6.666666667, 1.833333333, 2.833333333, 5.833333333, 5.166666667, 
2.941176471, 2.573529412, 6.25, 0.166666667, 0.666666667, 0.166666667, 
0, 0, 0.166666667, 0, 0, 0, 0.5, 0, 0.333333333, 0, 0, 0.333333333, 
0, 0, 0, 0, 0, 0, 0, 0.333333333, 0, 0, 0, 1.833333333, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0)

2个回答

您使用了错误的频率 $\omega$ , 或波长 $\lambda$ . 正弦函数的周期（“波长”）定义为 $\sin(\omega(t+\lambda)) = \sin(\omega t)$ ，因此

ω (t + λ) = 2 π ⟺ ω = \frac{2 π}{λ}

$\omega(t+\lambda) = 2\pi \iff \omega=\frac{2\pi}{\lambda}$ 您发布的数据中波长的粗略（视觉）猜测是 25。请注意，这不是您第一个图中显示的b $ zt的范围！如果我按如下方式拟合模型，我会得到一个不错的结果：

lmfit <- lm(walking ~ sin(2*pi/25*zt) + cos(2*pi/25*zt), data=b) 
plot(b$zt, b$walking)
x <- sort(b$zt)
lines(x, predict(lmfit, data.frame(zt=x)), col="red")

您可以使用最小二乘拟合从数据中自动估计周期（“波长”）。虽然 $\lambda$ 是一个非线性参数，因此不能直接由计算lm，您可以使用返回的残差lm来计算平方和并将其最小化 $\lambda$ ：

# function to be minimized
sum_squares <- function(lambda) {
  lmfit <- lm(walking ~ sin(2*pi/lambda*zt) + cos(2*pi/lambda*zt), data=b)
  return(sum(lmfit$residuals^2))
}
lambda <- optimize(sum_squares,c(1,50))$minimum

这计算 $\lambda\approx 26.05$ ，所以猜测 25 并没有那么糟糕。

您可能想做一些比拟合正弦波更灵活和更具探索性的事情，因为生物学不是物理学。

这是使用具有三个结的圆形样条拟合的广义加法模型 (GAM)（即，与正弦波的复杂性相当）。我们不应该使用最小二乘模型，因为它严重违反了它的假设：在每一天中，有很长一段时间几乎没有活动，而其他时间则有大量（倾斜）的活动分布。

要获得更详细的配合，请增加结数。 这是12节：

您可以通过这种方法了解更多关于这些数据的信息。即使您的最佳拟合碰巧看起来是正弦曲线，您也会找到证据证明它应该是正弦曲线。

我使用R包mgcv来创建这种配合，使用其内置的圆形样条（类型“cc”）。忽略它对“非整数 x”的看法：这个模型也适用于非计数数据。

library(mgcv)
# Assemble the data.  Notice the computation of `Time` as hour of the day.
X <- data.frame(zt = zt, Time=zt %% 24, Value=walking)

# Fit the model.  `k` is the number of knots.
fit <- gam(Value ~ s(Time, bs="cc", k=12), X, family="poisson", knots=list(Time=c(0,24)))

# Compute the predicted values for plotting.
Y <- data.frame(zt = seq(0, max(X$zt), length.out=201)[-1])
Y$Time <- Y$zt %% 24
Y$Prediction <- predict(fit, newdata=Y, type="response")

# Plot the data and the fitted values.
with(X, plot(zt, Value, main="Data with Detailed Fit")) 
with(Y, lines(zt, Prediction, lwd=2, col="Blue"))

其它你可能感兴趣的问题

上一篇当您拥有全部人口时，是否适合在数据上放置“误差线”？下一篇多峰分布的箱须图