2014-10-13

線形回帰に対する帰無仮説としてのウィーナー過程

一次線形回帰は $y~ax+b$ からの正規乱雑項で観測データを説明しようとしたときにa,bの最尤推定をすることであり、そのモデルの当てはまりのよさを、a=0という帰無仮説条件に照らして、検定することがよくあるし、aの推定値の信頼区間を考えるときも、a=0が信頼区間に入るのかどうかを気にするのが通例
それは、yがxによらずに一定であるところに乱雑項があるという条件で得られたデータであるという仮説を帰無仮説とする、ということ
回帰するのは、観測点の間に「つながりがない」ことを前提としているのだが、そんなことは気にせず、トレンドを見るために「線形回帰直線」をちょっと引いてみたりすることはよくある。そんな「回帰線」がどれくらい信用できるか、ということについての「めど」を与えるために…
帰無仮説がウィーナー過程である、ということにするのはどうだろうか？
ウィーナー過程でよく引き合いに出されるのが株価だが、株価のトレンドグラフを見て、傾向なんて存在していないのに、どうしても「傾向」を見てしまうのは、トレンドデータに対して、回帰直線(曲線)を想定してしまうこと、そのときに引き比べているのが水平線であること、水平線に比べれば、回帰直線(曲線)の尤度が十分に高いがために、「傾向あり」と感じるてしまうこと、と考えると、「過剰に傾向を読み取る」ことの説明がついて、本当の「傾向」は「ウィーナー過程」で説明をするよりも有意に回帰モデルの尤度が高いことを示すことによって確信するべきである、ということにはならないだろうか、ということ
さて。やってみる
まずはウィーナー過程でデータを作る。いろいろできるけれど、結構なトレンドのができたとする

# 観測時刻をばらつかせて作る
t <- sort(runif(n.pt))
diff.t <- diff(t)
k <- 1
# ウィーナー過程としてのx
x <- rep(0,length(t))
x[1] <- runif(1)
for(i in 2:length(t)){
	x[i] <- x[i-1] + rnorm(1,0,sqrt(k*diff.t[i-1]))
}
# x.wに格納しておく
# 一次線形＋正規乱雑項モデル
#x <- 3*t+2 + rnorm(length(t),0,0.1)
plot(t,x,type="l")
diff.x <- diff(x)

ウィーナー過程だとすると、単位時間当たりの分散を決めるパラメタを選ばないといけないので、それをグリッドで選ぶ

# ウィーナー過程の単位時間間隔を決めるパラメタをグリッドで作って、最尤推定値を求める
ks <- seq(from=0,to=k*2,length=100)
ks <- ks[-1]
ls <- ls.all <- rep(0,length(ks))
for(i in 1:length(ks)){
	ls[i] <- sum(-log(2*pi*ks[i]*diff.t)/2) - sum(diff.x^2/(2*ks[i]*diff.t))
	#ls.all[i] <- sum(-log(2*pi*ks[i]*diff.t.all)/2) - sum(diff.x.all^2/(2*ks[i]*diff.t.all))

}
plot(ks,ls)
abline(v=ks[which(ls==max(ls))],col=2)
max(ls)

最大尤度を出す時間単位のときの最大対数尤度

> max(ls)
[1] 2335.047

ついで一次線形回帰をして、その予測値と観測値との差の分散を推定し、その分散に相当する正規分布誤差が出ているものとして尤度を計算してみる

# 一次線形回帰で時刻tにおけるモデル予測値を出す
pr.lm <- predict(lm(x ~ t))
# 観測誤差の不偏分散
v. <- var(pr.lm-x)
# その分散を仮定して正規乱雑項の対数尤度
sum(-log(2*pi*v.)/2) - sum((x-pr.lm)^2/(2*v.))

(-log(2*pi*v.)/2) - sum((x-pr.lm)^2/(2*v.))
[1] -498.5286

あまりに「最大尤度」が違い過ぎて、不安になるが、ウィーナー過程仮説の対数尤度が大きく正であるのに対して、線形＋正規分布誤差モデルの対数尤度が大きく負であるから、「大小関係はＯＫ」の模様

x <- 3*t+2 + rnorm(length(t),0,0.5)
plot(t,x,type="l")
diff.x <- diff(x)
ks <- seq(from=0,to=k*2,length=100)
ks <- ks[-1]
ls <- ls.all <- rep(0,length(ks))
for(i in 1:length(ks)){
	ls[i] <- sum(-log(2*pi*ks[i]*diff.t)/2) - sum(diff.x^2/(2*ks[i]*diff.t))
	#ls.all[i] <- sum(-log(2*pi*ks[i]*diff.t.all)/2) - sum(diff.x.all^2/(2*ks[i]*diff.t.all))

}
plot(ks,ls)
max(ls)

では、線形＋正規乱数誤差のモデルが真である場合をやろう。今度は、線形＋正規誤差モデルの方の最大対数尤度に大きな違いはないが、ウィーナー過程モデルの方は、対数尤度の正負が逆転していて、ウィーナー過程モデルは、線形モデルに比べ大きく

> max(ls)
[1] -1150062

pr.lm <- predict(lm(x ~ t))
v. <- var(pr.lm-x)
sum(-log(2*pi*v.)/2) - sum((x-pr.lm)^2/(2*v.))

> sum(-log(2*pi*v.)/2) - sum((x-pr.lm)^2/(2*v.))
[1] -499.7419

図を二つ並べてみよう

上は、ウィーナー過程(酔歩)であって、下は直線的に右肩上がり、というのが「真実」
この２つの図をこのように別の過程の結果であると受け取るためには、どういう風にこの「絵」を読んでいればよいのだろう？