muttenz's blog

スイス星空だより

MCMC法でKaiV20の要素を改良 つづき

KatさんがMCMC法でEAのタイプの周期を求めることをすすめて下さり、この間ブログに書いたようにやってみました。

MCMC法でKaiV20の要素を改良 - muttenz's blog

そのスクリプトの

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])

この部分はEAのphase plotをいわば食外の水平な直線と、減光、極小、増光のV字の直線のモデルで表現しています。

つまりフェーズが0.0から前後p[3]の範囲内なら光度はp[4]+(p[3]-fph)*p[5]で表されるV字になり、それ以外の範囲では、食外の光度p[4]で一定光度になります。

Katさんが副極小もモデルに入れてみたらと書いておられたので、自己流にスクリプトの一部を修正してみた。もう一つV字をモデルに入れるわけです。 

yxep <- ifelse(fph < p[3], p[4]+(p[3]-fph)*p[5],

ifelse*1* p[7], p[4]))

「このブログでifelseのあとの部分にどうやってもp[4]+(p[6]-(0.5-fph))* p[7]が記入できず、へんな1になってしまう。どうしてだろう?」

2行目、p[4]+(p[6]-(0.5-fph))* p[7]が副極小のV字に相当し、最後のp[4]が食外光度になります。副極小のフェーズもパラメータとして入れることも可能とは思いますが、ここでは、簡単に副極小はフェーズ0.5で起こると言うモデルです。

これを組み込んだスクリプトで性懲りもなく、またKaiV20をやってみました。

パラメータをかなりきちんとやらないとうまく行きません。

得られたphase plot、見た目にすごく変わったわけではありませんが、微妙に数値が変わりました。

f:id:muttenz:20150125054643p:plain

下が前回、副極小なしのモデルでMCMC法を使って得たphase plot。

f:id:muttenz:20150115024953p:plain

今回、勉強していてパラメータの意味がはっきり分かりました。

paraminit <- c(P,E0,a,b,c)

Pは周期(日)、E0は元期-2450000、aは主極小の食の長さをd(日)とした時

a=(d/P)/2(食の長さをフェーズで表した値の1/2)、bは食外の光度m0、cは極小光度をm1とした時c=(m1-m0)/a

副極小を考慮した場合はパラメータが更に2つd、eと増え、dはaに、eは cに相当します。

勿論、これらの値を極めて正確に入れなくても計算を繰り返していく段階でこれらのパラメータがどんどん改良されていきます。しかし、とんでもない値を入れると時間はかかるしばあいによっては全く良い結果がでなくなる恐れがあります。なので、一応それほど外れていないと思われる値を入れたほうが時間の節約になります。aは勿論0.5以下。cはm1-m0>0、a>0なのでプラスです。

 

*1:0.5-fph) < p[6], p[4]+(p[6]-(0.5-fph