MCMC法でKaiV20の要素を改良
KatさんがvsoljのMLでEAの周期を正確に求めるには、サインカーブをモデルにしているLassoは向かず、MCMC法の方が良いとすすめてくださいました。
メーリングリストの記述を参考に、スクリプトをまとめました。(下の方に載せます。)
データはHJD、Magでdatに入れます。
最初のパラメータは今まで得られたデータを元に大雑把に入れます。
今回は paraminit <- c(1.7,7015.33,0.15,15.4,-2) としました。
1.7は大雑把な周期(day)、7015.33は元期から2450000を引いたもの、0.15は食の長さ(day)、15.4は食外の大雑把な光度、-2は食の減光の傾き(単位はmag/day)
paramsdには上記のパラメータの大雑把なsdを入れますが、ヤマカンで,
paramsd <- c(0.01,0.01,0.2,0.2,1) としましたが、これは計算を繰り返しているうちに改善されてきます。
これで入れなければならないパラメータは終わりで、これからスクリプトを最後まで動かします。
すると、plot(pp)でグラフが出ます。
そして第一回目のpとparamsdの値が表示されます。
> p <- apply(pres[,-(1:3000)],1,mean)
> p
[1] 1.64796743e+00 7.01533205e+03 8.24797264e-02 1.54462098e+01 2.23203816e+00
> paramsd <- apply(pres[,-(1:3000)],1,sd)
> paramsd
[1] 0.01603981941 0.00774792929 0.01476299806 0.01202531466 0.65402740705
それぞれ、最初に入れたパラメータやsdの値が変化しています。
第二回目の計算はスクリプトの途中、# after the first round, repeat from here のところから下を動かします。
またグラフやパラメータの値が表示されます。
2回めと同じように3回めもこの操作を続けます。3回めのグラフは
4回目のグラフは
このグラフがグチャグチャになったらほぼ目標達成ということだそうです。その時のパラメータの値は
> p <- apply(pres[,-(1:3000)],1,mean)
> p
[1] 1.74277168e+00 7.01533729e+03 5.35775292e-02 1.54276922e+01 7.01413228e+00
> paramsd <- apply(pres[,-(1:3000)],1,sd)
> paramsd
[1] 5.60280321e-05 2.28973080e-03 1.07440172e-03 1.79196337e-03 2.72777037e-01
でも念の為にもう一回やってみました。
そしてパラメータは
> p
[1] 1.74276935e+00 7.01533709e+03 5.33174946e-02 1.54281659e+01 7.07181988e+00
> paramsd <- apply(pres[,-(1:3000)],1,sd)
> paramsd
[1] 5.26242104e-05 2.08811145e-03 1.03913759e-03 1.81231895e-03 2.71882415e-01
各パラメータの値もsdももうあまり変化しないので、これで良さそうです。
ここで得られた元期2457015.33709と周期1.74276935 d,でphase plotを作ってみました。ちなみに以前PDMを使って得た周期は1.7427426 d.で元期はデータの中から言わばこの辺りという感じで選んでました。
見事完成!!それぞれのパラメータのsdが得られるのも素晴らしい。
Rのスクリプトを載せておきます。vsoljに載っていたのをところどころ順序を変えました。Katさん、こんなのでよろしいのでしょうか?
# vsolj 11049 2012.04.09
dat <- read.table("data.txt")
x <- dat$V1-2450000
y <- dat$V2
options(digits=9)
paraminit <- c(1.7,7015.33,0.15,15.4,-2) # (P,epoch-2450000,duration eclipse,max-mag,decline to min)
paramsd <- c(0.01,0.01,0.2,0.2,1) # sigma von ( , , , , )
magerr <- 0.1
nparam <- length(paraminit)
fac <- 1.0
nloop <- 10000
p <- paraminit
pres <- matrix(0,nrow=nparam,ncol=nloop)
pp <- numeric()
sq2 <- sqrt(2)
prob <- function(x,y,p) {
fph <- 0.5 - abs((((x-p[2])/p[1]) %% 1) - 0.5)
yxep <- ifelse(fph < p[3], p[4]+(p[3]-fph)*p[5], p[4])
return(-sum(((y-yxep)/(sq2*magerr))^2))
}
# after the first round, repeat from here
for (i in 1:nloop) {
pnew <- p + rnorm(nparam)*paramsd * fac
pr1 <- prob(x,y,p)
pr2 <- prob(x,y,pnew)
r <- pr2 - pr1
if (r > log(runif(1))) {
p <- pnew
pp[i] <- pr2
} else {
pp[i] <- pr1
}
pres[,i] <- p
}
plot(pp)
p <- apply(pres[,-(1:3000)],1,mean)
p
paramsd <- apply(pres[,-(1:3000)],1,sd)
paramsd