Transcript PPTX

かなた望遠鏡→
東広島天文台 外観↓
ベイズモデル
～最近の仕事から成分分離やトモグラフィーなど～
植村誠 （広島大学）
＠京大宇物談話会 110112
今日の話の流れ
• ベイズ統計とは
• ベイズモデルを使った最近の仕事例
– 矮新星のアウトバースト頻度を考慮した軌道周期分
布の推定
– ブレーザーの可視偏光データを利用した長期・短期
変動成分の分離
– 矮新星の早期スーパーハンプを利用した降着円盤
の高さマップの再構築
• まとめ
ベイズ統計とは
ベイズの定理 (Bayes’ theorem)
• データA,モデルパラメータBに対して、条件付確率
𝑃 𝐴 𝐵 𝑃(𝐵)
𝑃 𝐴 𝐵 𝑃(𝐵)
𝑃 𝐵𝐴 =
=
=
𝑃(𝐴)
𝑃
𝐴
𝐵
𝑃(𝐵)
𝐵
𝐿 𝐵 𝐴 𝑃(𝐵)
𝐵 𝑃 𝐴 𝐵 𝑃(𝐵)
– P(B): パラメータBの事前確率 (prior probability distribution)
– P(A|B) = L(B|A) : パラメータBの場合にデータAが生じる確率（＝尤度;
likelihood）
– P(B|A) : データAが観測された場合のパラメータBの事後確率
(posterior probability distribution)
– （Aが観測結果（定数）なのでP(A)は定数）
• 1763年、T. Bayesが発表。
• データから知りたいパラメータの「確率密度分布」が得られる。
– でも「事前分布」ってなんじゃ？
頻度主義とベイズ主義
頻度主義
ベイズ主義
真のモデルパラメータを推定する
モデルパラメータの事後確率分布を決定
する
データは不完全なサブサンプルの１つ
得られたデータが全て
データの数を増やすことで真の値に迫る
ことができる
データの数を増やすと事後分布の幅が
小さくなる
点推定＝最尤法（最小二乗法)
事後分布の平均or中央値
区間推定＝検定理論（χ2検定）から信頼
区間
事後分布の分散or中央値の前後で95%
が含まれる区間
事前分布ってなんやねん！怪しい！
事後分布の推定などできない！
（イメージです）
既知の情報を事前分布として利用できる。
無情報事前分布なら実際は最尤法と同
じ。
統計学者や科学哲学者ではない「１ユーザー」としては、「道具」としての有用性が
高くて理論がまともなら、主義は何でも良いでしょう。
でもベイズの定理の分母の積分項が解析的に解けない場合はどうする？
ベイズモデルと
マルコフ連鎖モンテカルロ法(MCMC)
• モンテカルロで事後分布の推定
– 計算機の進歩と共にここ数十年で急速な進歩
• MCMCの１つ、Metropolis法のアルゴリズム
– パラメータの組を x={x0,x1,….,xn} , 推定したい確率分布を P(x) として、
– xi から xi+1の候補を以下のように生成
𝒙𝑖+1 = 𝑁(𝒙𝑖 , 𝝈) : N(μ,σ)は平均μ,分散σ^2の正規分布
– P(xi+1)>P(xi) なら候補を採択
– P(xi+1)<P(xi) なら、確率 P(xi+1)/P(xi) で採択
– それ以外の場合は xi+1 = xi
尤度
最尤法で使われる
アルゴリズム。
間違ったピークを
捕まえることも
パラメータ
事後確率密度
MCMC。
分布を推定。多峰型
の分布にもメリット。
パラメータ
実例：滑らかな曲線モデル
（岩波書店 「階層ベイズモデルとその周辺」 石黒真木夫、他著 より）
•
•
課題：時系列データ a={ai} (i=1-10) をうまい
こと通るような滑らかな曲線 b={bj} (j=1-50)
を推定する。
モデル
– 𝑃 𝒃 𝒂 ∝ 𝑃(𝒂|𝒃)𝜋(𝒃)
– 𝑃 𝒂𝒃 =
1
𝑖 2𝜋𝜎
1
−(𝑎𝑖 −𝑏𝑖 )2
} （尤度関数）
2𝜎 2
−(𝑏𝑗 −2𝑏𝑗−1 +𝑏𝑗−2 )2
exp{
2
– 𝜋(𝒃) = 𝑗
exp{
}
2𝜔2
2𝜋𝜔2
（事前分布：２階階差数列が正規分布に従う）
宇宙物理学におけるベイズモデルの応用
論文数
（ADSでアブストラクトに”Bayes”の文字列が含まれる査読論文）
300
250
200
150
100
50
0
1980
1985
1990
1995
2000
この10年間で急速に普及
2005
2010
ベイズモデルの利点・弱点
•
利点
– 既知情報を事前分布としてモデルに組み込める
– MCMCとの相性の良さ
•
•
逆問題も解ける
パラメータスペース総当たりのモンテカルロよりは「次元の呪い」に強い
•
MEMなど従来の手法よりも自由度が高い
–
高次元の画像再構成など。
– 宇宙物理に向いている面も
•
•
何度も繰り返し「実験」ができない、光子が少ない、観測機会が少ない、情報が縮退してる、などなど。
弱点（注意点）
– 適切な尤度関数（モデル）の設定が必要なことは通常の統計理論と同じ
– 事前分布の設定の妥当性・正当性
– 事前分布にハイパーパラメータが含まれる場合の扱い
•
階層ベイズ or 経験ベイズ
– MCMCの収束判定、サンプリングにもパラメータが必要
•
最近、ベイズ的な仕事をいくつかやってきたので、今日はそれらを紹介します。
ベイズモデルを使った最近の仕事例 その１
矮新星のアウトバースト頻度を考慮
した軌道周期分布の推定
(Uemura, et al., 2010, PASJ, 62, 613)
激変星とは
伴星（主系列星）
降着円盤
白色矮星
• 低質量連星系の１つの末路
軌道周期分布の問題
•
激変星の進化＝連星系からの角運動
量の抜き取り
– Magnetic braking
– Gravitational wave
•
最小周期
– 伴星のコアが縮退（brown dwarf）しはじめ
ると、星の「質量ー半径」の関係が変化
– 通常の伴星：短周期の系へ進化
– 縮退した伴星：長周期の系へ進化
– 観測される最小周期 ～８０分
•
連星の初期状態を仮定して population
synthesis
– 大半の激変星は最小周期付近にいるは
ず（period spike）。
– 最小周期（period minimum）は60－70分
•
観測される分布を説明できない
WZ Sge型矮新星
• 矮新星
– 降着円盤の熱的不安定
性で増光
– 増光間隔は数日～数年
2003年までの軌道周期分布と、2008年までの軌道周
期分布。最小周期付近で特に増加している。
• WZ Sge型矮新星
– 10年以上の増光間隔
– 軌道周期が最小周期付
近
– 最近のAll sky monitorな
どで発見頻度が増加
• 軌道周期分布が変化してきた →矮新星のアウ
トバースト頻度を考慮することが重要
観測された軌道周期分布から本来の分布をベイズ的に推定
•
軌道周期と増光頻度の関係
増光が観測される天体の軌道周期の確
率密度関数
本来の軌道周期分布を表わす
確率密度関数のモデル
爆発頻度に依存した発見確率
絶対等級に依存した発見確率
• 観測されるアウトバースト
＝ある確率分布からの標本、と考える
サンプル
• 2003年から2007年の5年間で見つかった矮新星アウトバースト
– ASASで42天体 （単一の観測システム）
– VSNETで146天体 （様々な観測システムが混合）
パラメータの事後分布
• 最小周期は観測値よりも短くなる
矮新星の軌道周期分布のベイズ的推定
まとめ
• 最近の観測から、WZ Sge型矮新星が従来よりも高い
頻度で発見。SDSSによる暗い（伴星からの質量輸送
率が低い＝最小周期付近の）激変星の検出も合わせ
て、軌道周期分布が変わってきた。
• 矮新星のアウトバースト頻度を簡単な関数でモデル
化して、ベイズ的に本来の軌道周期分布を推定。
• Period minimum問題、Period spike問題は WZ Sge型
星（＝増光頻度の低い矮新星）を考慮すれば解決で
きるかもしれない。
ベイズモデルを使った最近の仕事例 その２
ブレーザーの可視偏光データを利用
した長期・短期変動成分の分離
(Uemura, et al., 2010, PASJ, 62, 69)
blazar
ブレーザーとは
Optical-NIR
High polarization
A probe of the magnetic field structure
But, no universal law has been established
• ジェットの研究に貴重な天体
ブレーザーの偏光観測
Positive correlation between the
flux and PD. (Smith, et al. 1986)
Apparently random motion in the
QU plane (Moore, et al. 1982)
• 偏光の変動に規則性はあるのか？
偏光の２成分モデル
Stokes parameters for
linear polarization
• 長期変動と短期変動
に分離したい
Stokes’ QU parameters of BL Lac for ~20 yr
(Hagen-Thorn, et al. 2002)
(0,0)
偏光の成分分離のためのベイズモデル
U
p ( Q 0 , U 0 | f , Q obs , U obs ) 
Q
L ( f , Q obs , U obs | Q 0 , U 0 )   ( Q 0 , U 0 )
Posterior distribution
of the long-term trend
flux
t
Likelihood function,
maximized when the
total flux and polarized
flux completely
correlate.
C
Prior distribution of the
long-term trend
(smoother curve is
preferred).
 (Q 0 ) 

 ( Q 0 , i  Q 0 , i 1 ) 2 
exp  

2
2
2w
2 w


1
• 短期変動の偏光フラックスが総フラックスと相関
するように長期変動を推定する
人工データを使ったテスト
pol. flux
pol. flux
total flux
total flux
corrected pol. flux
long-term trend
(0,0)
(0,0)
U
U
• 期待通りの結果
Q
Q
人工データを使ったテスト
pol. flux
pol. flux
total flux
total flux
corrected pol. flux
long-term trend
(0,0)
(0,0)
U
U
• 期待通りの結果
Q
Q
人工データを使ったテスト
pol. flux
pol. flux
total flux
total flux
corrected pol. flux
long-term trend
(0,0)
(0,0)
U
U
• 期待通りの結果
Q
Q
人工データを使ったテスト
pol. flux
pol. flux
total flux
total flux
corrected pol. flux
long-term trend
(0,0)
(0,0)
U
U
期待通りの結果
• 期待通りの結果
Q
Q
実際の観測データでは
• OJ 287
• S2 0109+224
pol. flux
total flux
total flux
pol. flux
(0,0)
U
U
(0,0)
Q
Q
実際の観測データでは
• OJ 287
• S2 0109+224
pol. flux
total flux
total flux
pol. flux
長期成分の偏光方位角が振動？
長期的なジェット軸の運動か？
(0,0)
U
U
(0,0)
Q
Q
実際の観測データでは
• OJ 287
• S2 0109+224
pol. flux
total flux
total flux
pol. flux
長期成分の偏光方位角が振動？
長期的なジェット軸の運動か？
(0,0)
U
U
(0,0)
Q
Q
実際の観測データでは
• OJ 287
• S2 0109+224
pol. flux
total flux
total flux
pol. flux
長期成分の消長が見えている？
長期成分の偏光方位角が振動？
長期的なジェット軸の運動か？
(0,0)
U
U
どちらも興味深い結果に。
•
•
(0,0)
１シーズンだけでなく、複数シーズンの結果を見たい。
長期成分の挙動と電波の挙動とに相関があると面白いかも。
Q
Q
ベイズモデルを使ったブレーザー偏光の成分分離
まとめ
• ブレーザーの偏光に法則性が見られないのは
長期変動成分が短期フレアに重なっているから
かもしれない。
• その仮説に基づいてベイズ的に偏光成分を分離
するモデルを開発。
• いくつかの天体では２成分モデルで上手く説明
できることが可能で、「長期偏光成分の挙動」と
いう新しいアプローチを得た。
ベイズモデルを使った最近の仕事例 その３
矮新星の早期スーパーハンプを利用
した降着円盤の高さマップの再構築
(Uemura, et al., まだ解析中)
早期スーパーハンプとは
上：早期スーパーハンプと通常のスー
パーハンプの光度曲線（と色変化）。デー
タはV455 Andのもの。
下：WZ Sgeのアウトバースト光度曲線
• 円盤への潮汐効果 →ただし詳細不明
円盤の高さを見ている？
Early superhumps in V455 And (Matsui, et al. 2009)
ハンプ成分が赤い →円盤外縁が膨張？
Zoo of early superhumps (Kato , et al. 2002)
Edge-onの天体ほど振幅が大きい？
Two-armed spirals on the disk and
early superhumps
(Maehara, et al. 2007)
• 回転効果 →観測から円盤を再構成できるかも
早期スーパーハンプのモデル
•
The 2:1 resonance
–
Osaki & Meyer (2002)
Tidal truncation and resonance radii
•
Tidally distorted disk
–
–
Kato (2002)
Even below the 2:1 resonance radius
Height distortion in the disk (Ogilvie, 2002)
伴星
• 観測から円盤を再構成して理論モデルと比較
ベイズを使った円盤高さの再構成
𝑃 ℎ 𝑑 ∝ 𝐿 ℎ 𝑑 𝜋𝑠𝑚𝑜𝑜𝑡ℎ 𝜋𝑠𝑡𝑑.𝑑𝑖𝑠𝑘
(𝑓𝑖,𝑑𝑎𝑡𝑎 −𝑓𝑖,𝑚𝑜𝑑𝑒𝑙 )2
• 𝐿 ℎ 𝑑 ∝ 𝑖 exp(−
)
2𝜎 2
→モデルと観測の光度曲線で尤度
• 𝜋𝑠𝑚𝑜𝑡ℎ ∝
𝑖,𝑗 exp
−
(ℎ𝑖𝑗 −2ℎ𝑖−1,𝑗 +ℎ𝑖−2,𝑗 )2
2𝑤1
2
exp(−
(ℎ𝑖𝑗 −2ℎ𝑖,𝑗−1 +ℎ𝑖,𝑗−2 )2
2𝑤1
2
)
→隣り合うグリッドの高さはなるべく滑らかに（事前分布）
• 𝜋𝑠𝑡𝑑.𝑑𝑖𝑠𝑘 ∝
𝑖,𝑗 exp(−
(ℎ𝑖𝑗 −ℎ𝑠𝑡𝑑.𝑑𝑖𝑠𝑘 )2
2𝑤2
2
)
→標準円盤（h=0.1r）を仮定（事前分布）
•
•
使う観測結果は、g,V,Rc,Ic,Jバンドのphase-averaged light curve (20 points/phase)
座標(r,θ)での高さhを推定する。
– 新しい h 候補を発生
– モデル光度曲線 →尤度計算
– Metropolis または Multiple-Try Metropolis 法で採択・棄却を判定
•
前２例よりも複雑なモデル、多い次元
人工データを使ったデモ
観測
再構成
• 外側内側が区別できている
矮新星 V455 And のデータを使った結果
(preliminary)
伴星
• アウトバースト極大から5日目の
データ
• 結果が何を語っているかはこれ
から検討
ベイズモデルによる降着円盤の再構成
まとめ
• 早期スーパーハンプは回転してる円盤の高さ
方向の非軸対称性を見ている。
• 多バンド測光観測から円盤の構造をベイズ
的に再構成するモデルを開発（中）
• 使い物になりそう。今年中に論文化予定。
– 早期スーパーハンプの多バンド光度曲線から、
円盤構造の時間変化が追える
今日の話のまとめ
• ベイズ統計＋MCMCを使った手法は宇宙物理の分野
でもここ10年で急速に広がった。
– 成分分離、逆問題、少ない観測点、などの状況で便利
• ３例を紹介
– 矮新星の本来の軌道周期分布の推定
– ブレーザーの偏光成分の分離
– 早期スーパーハンプから円盤構造の再構成
• 1000次元以下なら収束する実績あり。アイデアのある
方、いっしょに勉強しませんか？
ありがとうございました。
予備トラぺ
伴星
推定された円盤。極座標でグリッドを設定して、かつ動
径方向は対数スケールで等間隔にしてる。普通のソフト
だと表示しにくいので線形スケールで等間隔に補間した。
WZ Sgeアウトバースト初期のHeIIのドップラー
トモグラフィー(Baba, et al. 2002)
今回のデフォルト
外縁温度を1割高く
温度勾配を標準円盤よりも緩く
（p=-0.75→-0.65）
Inclination = 80 deg.
等級
ここのギャップは不明。赤側では目立たないので、定
性的には円盤の比較的外側がかなり膨らんでいて、内
側の高温領域を隠せば良い。ただしあまり膨らみ過ぎ
るとそもそも暗くならない。
光度曲線の比較。赤が観測
で、青が再構成されたもの。
図中のphase=0.2くらいの極小
が、orbital phase = 0.0（伴星が
手前）。
観測値の誤差は一律0.01を仮
定した。
ここのギャップは伴星による円盤外縁部の食
で説明できそう。Inclination~75degだと、内側
までは隠せない。
横軸：ＭＣＭＣのステップ、縦軸：事後確率密度。
混合がよくない。
ＭＣＭＣサンプリングの違い。
左：3万点から10万点まで100点ずつ。右：3万点から100万点まで100点ずつ。
スケールもあるが、特に内側の詳細構造が見た目に違う？
点数が多い右の方がよりスムーズ。ただし、正しい結果を出すためにはまず混合を
改善しないといけない。
設定したパラメータの詳細
• V455 And のパラメータ (Araujo-Betancor, et al. 2005)より
– Incliantion = 75deg
• 浅い食があるため。精度は良くないはず。
– Mass ratio~9 (M1~0.6Mo, M2~0.07Mo)
– Porb = 0.056309 d
• モデルに使ったパラメータ
– Inclination = 75 deg
– Rin=1(WD半径）として、Rout=30（2:1 resonance に対応）
– 外縁温度 6800K
• 観測された日の平均の色から推定
– 温度は T∝r^(-0.75) （標準円盤）、ただし、r = sqrt(x^2+y^2)、つまり高さの寄
与は無視。
• 計算速度を上げるため。
– 動径方向16分割、方位方向20分割
– Irradiationや伴星の食の効果は入っていない