時系列データのベイズ推定アルゴリズム(モクモク開発)

とりあえずMemo：

\displaylines{
Z(t) \sim N(HX(t),{\sigma}_{1}^{2}I), X(t)\sim N(FX(t-1),I)を考えるとき, \\
単純のため, X(t)=\sum_{k=1}^{m} a_{k} X(t-k)+w, w \sim N(0,1)\\
           Z(t)=\beta X(t) +v, v\sim N(0,{\sigma}_{1}^{2}) \\
として, 各パラメータが\beta \sim N(1,1), (a_{1}, ... , a_{m}) \sim N(0,{\sigma}_{1}^{2}I),\\
{\sigma}_{1} \sim Invgamma(\sigma, 3, 2) 
に従うとして, 計算すると, 尤度部分が\\
\frac{N}{\sqrt{2 \pi (\beta^{2}+\sigma_{1}^{2})}} \prod_{t}^{N} exp(-\frac{(Z(t)-\beta \sum_{k=1}^{m} a_{k} X(t-k))^{2}}{2(\beta^{2}+\sigma_{1}^{2})})\\
となる(はず). 各パラメータの分布から事前分布を考えてメトロポリスヘイスティ\\
ング法を行って計算した. Z(t), X(t)はデータとして得られている構想.\\
}

おっ、なんかとりあえず使えそう（笑）

回帰パラメータの事前分布の平均が０にしてあるから、尤度の影響が少ないパラメータは事前分布を優先して勝手に０になっているところがよろし。

data("EuStockMarkets")
dat=EuStockMarkets

Z=dat[,1]
X=dat[,2]

Y=as.numeric(Z)

#特徴量と目的変数を編集
p=6
y=c()
Xn=array(0,dim=c(length(Y)-p,p))
for(j in 1:nrow(Xn)){
val=Z[length(Y)-j+1]  
y=c(y,val)
vec=X[(length(Y)-j):(length(Y)-j-p+1)]
Xn[j,]=vec  
}
N=length(y)
Y=y



#確率勾配法にて最尤推定した値を初期値としてベイズ推定
#メトロポリスヘイスティング法
library(MCMCpack)
sigma=1

#逆ガンマ分布のパラメータ
v=3 #Scale
w=2*sigma #Shape
n=ncol(Xn)

posterior=function(x){
aa=x[1:(n)];x=x[-c(1:n)]
bb=x[1];x=x[-c(1)]
ss=x
likelihood=-length(Y)*log(2*pi*(bb^2+ss^2))/2-sum((y-bb*Xn%*%aa)^2)/(2*(bb^2+ss^2))
return(likelihood+sum(log(dnorm(aa,0,1)))+sum(log(dnorm(bb,1,1)))+sum(log(dinvgamma(ss^2,v,w))))
}


ite=6000
eta=10^(-1)
posterior_values=c()
A=(solve(t(Xn)%*%Xn)%*%t(Xn)%*%Y)
B=1
S=1

for(l in 1:ite){


for(j in 1:length(A)){
  
#Aの更新
a0=rnorm(1,A[j],1)  
A0=A;A0[j]=a0
r=min(c(1,exp((posterior(c(A0,B,S))-log(dnorm(a0,A[j],1)))-(posterior(c(A,B,S))-log(dnorm(A[j],a0,1))))))
val=sample(c(0:1),1,prob=c(r,1-r))
#if(val==0){A=A0}
#if(val==0){A=A*(1-r)+A0*r}
if(val==0){A=A*(1-eta)+A0*eta}
A[A<0]=10^(-5)
}

#Bの更新
b0=rnorm(1,B,1)  
B0=B;B0=b0
r=min(c(1,exp((posterior(c(A,B0,S))-log(dnorm(b0,B,1)))-(posterior(c(A,B,S))-log(dnorm(B,b0,1))))))
val=sample(c(0:1),1,prob=c(r,1-r))
#if(val==0){B=B0}
#if(val==0){B=B*(1-r)+B0*r}
if(val==0){B=B*(1-eta)+B0*eta}  
  
#Sの更新
s0=rnorm(1,S,1)  
S0=S;S0=s0
r=min(c(1,exp((posterior(c(A,B,S0))-log(dnorm(s0,S,1)))-(posterior(c(A,B,S))-log(dnorm(S,s0,1))))))
val=sample(c(0:1),1,prob=c(r,1-r))
#if(val==0){S=S0}
#if(val==0){S=S*(1-r)+S0*r}
if(val==0){S=S*(1-eta)+S0*eta}


print(posterior(c(A,B,S)))
posterior_values=c(posterior_values,posterior(c(A,B,S)))

}

plot(posterior_values)
cor(Y,B*Xn%*%A)