MCMC法でEAの要素を求める際の、パラメータやその入れ方を工夫してみた。
Katさんのご意見では、MCMC法でやる場合にパラメータの数が増えるとそれだけ正しそうな解になかなか到達しないので、EAの場合、副極小はフェーズ0.5と固定してもよいのではということです。
今までのスクリプトでは食外から極小への傾きがパラメータとして入っていましたが、それは大雑把でも毎回計算しないといけないので、もっと把握しやすいパラメータ、極小光度をパラメータにしました。
paraminitで数字が並ぶとどれがどれだかよくわからなくなるので、パラメータを入れる部分を前に置きました。
これらを工夫して書きなおしたスクリプトを下に載せておきます。
またKaiV20で試した結果。順に、
周期、 元期、 食外光度、 第一食の幅の半分、
第一極小光度、 第二食の幅の半分、第二極小光度。
1.74277642e+00 7.01533708e+03 1.54228110e+01 5.43491780e-02 1.58023035e+01 6.37331798e-02 1.55301130e+01
これだと、食外光度が15.42、第一食の幅が10.9％(5.43*2)で極小光度が15.80、第二食の幅が12.7％(6.37*2)で極小光度が15.53と必要な要素がすぐに見て取れる。
上記の要素のｓｄ。
2.95419457e-05 1.18810090e-03 9.01790171e-04 5.14007862e-04
4.41436032e-03 2.55451990e-03 5.88479472e-03
新しいスクリプト
# vsolj 11049 2012.04.09
# MCMC methode with min II parameter p[6],p[7]
dat <- read.table("data.txt")
x <- dat$V1-2450000
y <- dat$V2
# パラメータを入れる。sdはそれぞれのパラメータのsdの値（適当だが、
# 毎回ほぼ同じに出来る）
p <- 1.74 #period
sdp <- 0.05
epo <- 7015.33 # epoch - 2450000
sdepo <- 0.05
mag <- 15.43 # mag. out of ecliipse
sdmag <- 0.05
d1 <- 0.05 # the half of phase duration min1
sdd1 <- 0.03
min1 <- 15.8 # mag. of min1
sdmin1 <- 0.05
d2 <- 0.05 # the half of phase duration min2
sdd2 <- 0.03
min2 <- 15.5 # mag. of min2
sdmin2 <- 0.05
magerr <- 0.05
options(digits=9)
paraminit <- c(p,epo,mag,d1,min1,d2,min2)
paramsd <- c(sdp,sdepo,sdmag,sdd1,sdmin1,sdd2,sdmin2)
nparam <- length(paraminit)
fac <- 0.5
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[4], p[5]-(p[5]-p[3])*fph/p[4],
ifelse ( (0.5-fph) <=p[6], p[7]-(p[7]-p[3])*(0.5-fph)/p[6], p[3]) )
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
