Tokyo¶このノートでは、JpGU2018にて発表したTRAPPIST-1の惑星質量推定のためのライトカーブ解析についての方法をまとめました。この研究は2017年度卒業研究として行ったものです。測光には福井暁彦さんのプログラムafapphotを用いた。様々なアパーチャー半径で測光した結果、目で見てライトカーブが綺麗に見えるアパーチャー半径を複数用意した。以下で行うフィッティングを様々な半径で行なった結果BIC等の基準から最も良くフィッティングできる半径を残すことにする(結果としては13pixだった)。0.005028GJDをBJDに変換したものに置き換える。[5]:len(df)Out[5]:306データの数は306点。これをプロットしてみると以下のようになる。上のグラフから、1箇所明らかに外れている値があるのがわかる。確認したところすべておなじフレーム(MSCT2_1711050228)のデータだったので、fluxの平均から$5\sigma$以上離れている点を除くと、その1点だけが除かれた。5df=df[np.abs(df['flux(12)']-df['flux(12)'].mean())<=(sigma*df['flux(12)'].std())]また、airmassについても1点だけおかしい点があった。header書き込みエラーだと思われるので、データ点を減らさないために、前後の平均値を用いることにして使うデータに戻してやった。[10]:len(df)Out[10]:305外れ値を除いたデータ点は305点、プロットすると以下のようになる。0x1153bf080>左側の凹みがe、右側がdのトランジットと予想されるが、確認はモデルフィットの際に行う。)を、2つのトランジットが重なった状況にも適応できるようにした。フィッティングの際に半径比kとトランジット中心tcのみを変動させ、その他のパラメタは固定している。状況に応じて固定したパラメタをフリーにしたり、priorの置き方も工夫する。必要な物理量の代入¶トランジットモデルを描くための物理量を先行文献から求める。)を用いて周縁減光パラメタの理論値を求める。u2_,u2err周縁減光のパラメタとして、$(u_1,u_2)$よりも以下で定義される$(q_1,q_2)$の方が使いやすい(範囲を指定すれば物理的にありえない状況を除ける)ため、変換の式。u_to_q(ldc[0],ldc[2])ためしに生データに理論値を入れたトランジットモデルを重ねると以下のようになる。left')pl.tight_layout()本来はTRAPPIST-1dとeのトランジットが重なる予定だった(それを狙った観測であった)。実際にはTTVの影響でタイミングが重ならなかったことがわかる。$という線形の組み合わせが成り立っているとする。($X_{im}$はairmassやdxのような要素、w_mはパラメタとなる重み。pythonの行列計算を用いる。逆行列を求めるため正方行列にする必要があるので以下のようにしている。systematicsの組み合わせを用意する¶systematicsが上記のような線形の組み合わせで表されると近似した場合、すべての要素を入れれば入れるほど正しい結果になるとは限らず、どの要素を入れるかを考慮する必要がある。プログラムを用いてその組み合わせを用意する。sys_combination.append(X)このうちどの組み合わせを$X$として採用するかをフィッティングして(BIC等を基準にして)決めたい。$w$がフリーパラメタとなるが、初期値が適切でないとおかしな結果を出してしまうので、初期値として、OOT(Out-Of-Transit)からざっと求めた値を用いることにする。各systematicsの組み合わせについて$\vec{w}$の値を求める¶まず、OOTのデータ点(ざっと)を用意する。w_before.append(w)これで各systematicsの組み合わせについて$\vec{w}$のおよその値が求められた。各パラメタセットでフィットし最も良いセットを決める(MLE)¶systematicsについて様々なパラメタの組み合わせを用意した。パラメタは、増やせば増やすほどよくフィットされてしまうが、これは物理的に正しくない「overfitting」を起こしている可能性がある。そこで、カイ自乗でなくBICという値を用意し、これを小さくするようにフィッティングを行う。最小にする関数を定義。np.sum(((f-s*m)/err)**2)使う(最適化する)パラメタを整える。簡単にフィッティングしてみるため、numpy.optimize.minimizeをつかう。298.1162optimize後のモデルは以下のようになった。model');pl.legend(fontsize=18)fig.tight_layout()さらにここで、(短いチェーンの)MCMCを用いて、optimizeしたものからさらに少し動いた位置での値を試してより最適な方を選択する、という手法を導入する。いうパッケージを用いた。294.8089MCMCを使った場合の方がBICが3小さくなった。model');pl.legend(fontsize=18)fig.tight_layout()文献値、optimizeの結果、mcmcによるフィッティングの結果を並べると以下の通り。-4.03808313e-03]])同様のことを様々なアパーチャー、31パターンのパラメタの組み合わせで行い、BICを計算する。(radiusの値を変化させる)なお、errはreduced$\chi^2$がおよそ1になるように0.935倍した(測光プログラムでおよそ良いエラーを出しているが過大・過少見積もりしている場合があるのでここで補正する)。348.5125ここで様々な半径を試したところ、13pixが最も良いという結果を得た。これで全てのパラメタセットの場合のBICが計算できた。これを小さい順に並び替える。0.979682このあとMCMCでつかいたい、BICが最小の10組のパラメタセットをまとめておく。良いパラメタセットでフィッティングする(MCMC)¶emceeについてはこれをみるべし(https://speakerdeck.com/dfm/data-analysis-with-mcmc)以下のnumを0~9まで変化させるとそれぞれのパラメタセットについてのMCMCを実行できる。[49]:free[good_number[num]]Out[49]:['airmass']nstepsを変えることでより精度を上げられる。lengthを計算し、それによりステップ数を調節して、収束するのに十分なステップ数を与えている。326.2667BICもわずかによくなっている。これを答えとする前に、本当にMCMCのチェーンは収束したのかを確認する必要がある。いくつかの収束の確認法を行う。各ステップごとのチェーンの動きを見てみる。(100点が10000ステップ動いているような図。0x118132da0>チェーンを眺めてみると、程よく広がったところから収束していっているのを確認できる(そうなるように初期値のばらつきを設定した。)チェーンの最初の部分はまだ収束していないので(burn-inと呼ばれる)、ここは取り除いて各パラメタの分布を見てみる。fontsize=10)。