EAの副極小のフェーズもパラメータにしてみた

下のスクリプトのようにやってみた。p[8]は副極小が0.5からどれだけずれているかというパラメータ。

fphの代わりになんとかph <- {(x-p[2])/p[1]}%%1だけで出来ないかと

abs{ph-round(ph)}とやってみた。しかしこれで計算させたらどうしても数値がおかしくなる。どうもコンピュータがやる数値の丸め方が数学的にroundするのとは異なるらしい。

そこで、下のスクリプトではfphと phが混在してしまった。

しかし、どうもうまく働かない。特に周期が外れていってしまう。原因が良くわからない。paramsdをあまりおおらかに取るとうまくいかない感じがする。カッコはRでの使い方のように書くとこのブログでは変になるのでところどころRの書式からは外れています。

# vsolj 11049 2012.04.09
# with min II parameter p[6],p[7],p[8]

dat <- read.table("data.txt")
x <- dat$V1-2450000
y <- dat$V2

options(digits=9)
paraminit <- c(1.74, 7015.33, 0.05, 15.4, 7,0.05,2,0)
# (P,epoch-2450000,duration minI/2,max-mag,decline to minI,
# duration minII/2,decline minII,dif. from minII)
paramsd <- c(0.005, 0.01, 0.01, 0.07, 1,0.02,0.5,0.001) # sigma von ( , , , , )
magerr <- 0.1

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) {
ph <- {(x-p[2])/p[1]}%%1
fph <- 0.5 - abs(ph - 0.5)
yxep <- ifelse(fph < p[3], p[4]+(p[3]-fph)*p[5],
ifelse(abs(ph-0.5-p[8])< p[6], p[4]+(p[6]-abs(ph-0.5-p[8]))* p[7], 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