まあ，
ぼちぼち考えてみましょう．
上の試行錯誤でつかった R スクリプトはこんなかんじで．
n.rep <- 30
n.year <- 50
sigma.m <- 0.1
sigma.n <- 0.8
mean <- replicate(n.rep, cumsum(rnorm(n.year, 0, sigma.m)))
noise <- replicate(n.rep, rnorm(n.year, 0, sigma.n))
y <- mean + noise
ylim <- c(-1, 1) * (sigma.m * sqrt(n.year) + sigma.n) * 2.5
par(mfrow = c(2, 1))
par(mar = c(2.5, 2.5, 0.5, 0.5), mgp = c(1.5, 0.5, 0))
matplot(mean, type = "l", xlab = "x", ylim = ylim)
matplot(y, type = "l", xlab = "x", ylim = ylim)
fit.glm <- function(yy, xx)
{
m1 <- glm(yy ~ 1)
m2 <- glm(yy ~ xx)
ifelse(m1$aic < m2$aic, 1, 2)
}
xx <- 1:n.year
print(table(sapply(1:n.rep, function(i) fit.glm(y[,i], xx))))
