はじめに
以前、X線天文学で固有運動を測るときに使われるプロファイルシフトの統計量について、以下の記事を書きました。
そこで紹介したのが、Katsuda et al. (2008) などで用いられている次の形です。
Q(s)
=
\sum_i
\frac{
[f_i-g_i(s)]^2
}{
\sigma_{f,i}^2+\sigma_{g,i}^2(s)
}.
ここで $f_i$ と $g_i(s)$ は異なる時期に得られたプロファイル、$s$ は試行する位置ずれです。
以前の記事では、
「観測側にも比較するプロファイル側にも誤差があるので、分散を足して考える」
という直感的な説明をしていました。
これは二乗項を理解するには分かりやすいのですが、最近もう少し統計的な出発点から考え直してみると、
そもそも $f$ も $g$ も観測値なのだから、最初から2つの観測の同時尤度として書く方が自然ではないか？
ということが気になりました。
さらに考えている途中で、
プロファイルをずらすと比較範囲に入るbinが変わる。
ではGaussianの正規化項も位置ずれ $s$ に依存するのではないか？
という疑問も出てきました。
今回は、このあたりを整理したメモです。
1. 普通のGaussianカイ二乗
まず通常の問題を考えます。
観測値 $y_i$ が
y_i
\sim
\mathcal{N}
\left(
\mu_i(\theta),
v_i
\right)
に従うとします。
Gaussian尤度は
L(\theta)
=
\prod_i
\frac{1}{\sqrt{2\pi v_i}}
\exp
\left[
-\frac{
[y_i-\mu_i(\theta)]^2
}{
2v_i
}
\right].
したがって、
-2\ln L(\theta)
=
\sum_i
\frac{
[y_i-\mu_i(\theta)]^2
}{
v_i
}
+
\sum_i\ln(2\pi v_i).
$v_i$ が推定したいパラメータ $\theta$ に依存しなければ、第2項は定数なので、
\chi^2(\theta)
=
\sum_i
\frac{
[y_i-\mu_i(\theta)]^2
}{
v_i
}
を最小化すればよいことになります。
普通の「観測値とモデルを比較するカイ二乗」はこの形です。
2. 固有運動では2つとも観測値
一方、固有運動のプロファイル比較では少し状況が違います。
2つの観測プロファイルを
f_i,
\qquad
g_j
とします。
それぞれの分散を
a_i=\sigma_{f,i}^2,
\qquad
b_j=\sigma_{g,j}^2
とします。
ある位置ずれ $s$ を試したとき、$f_i$ と対応する $g$ のbinを $g_{j_s(i)}$ と書きます。
ここで、一方を「観測値」、もう一方を「誤差を持つモデル」と考える必要はありません。
むしろ、
$f_i$ と $g_{j_s(i)}$ は、同じ未知の真の強度 $m_i$ を独立に観測したもの
と考えることができます。
つまり、
f_i
\sim
\mathcal{N}(m_i,a_i),
g_{j_s(i)}
\sim
\mathcal{N}
\left(
m_i,
b_{j_s(i)}
\right).
です。
これを2つの観測に対する同時尤度の出発点とします。
3. 共通の強度 m_i を最適化する
1組の対応するbinについて、$-2\ln L$ の二乗項は
\frac{(f_i-m_i)^2}{a_i}
+
\frac{
[g_{j_s(i)}-m_i]^2
}{
b_{j_s(i)}
}.
知りたいのは $m_i$ ではなく位置ずれ $s$ なので、$m_i$ は補助的な未知量として最適化します。
最適値は
\widehat m_i
=
\frac{
f_i/a_i+
g_{j_s(i)}/b_{j_s(i)}
}{
1/a_i+
1/b_{j_s(i)}
}.
つまり、2つの観測値の逆分散重み付き平均です。
これを戻すと、
\min_{m_i}
\left[
\frac{(f_i-m_i)^2}{a_i}
+
\frac{
[g_{j_s(i)}-m_i]^2
}{
b_{j_s(i)}
}
\right]
=
\frac{
[f_i-g_{j_s(i)}]^2
}{
a_i+b_{j_s(i)}
}.
したがって、
Q(s)
=
\sum_i
\frac{
[f_i-g_{j_s(i)}]^2
}{
a_i+b_{j_s(i)}
}
が得られます。
これはまさにKatsuda型の統計量です。
つまり、
\frac{(f-g)^2}
{\sigma_f^2+\sigma_g^2}
という形は、
「両方に誤差があるから分散を足す」
という直感だけでなく、
2つの観測が共通の未知強度を測っていると考え、その未知強度を同時尤度から最適化して消去した結果
としても導くことができます。
4. 残差から考えても同じ二乗項は出る
ここで別の考え方もあります。
残差を
r_i(s)
=
f_i-g_{j_s(i)}
と定義します。
$f_i$ と $g_{j_s(i)}$ が独立なGaussianなら、固定された1組については
r_i
\sim
\mathcal{N}
\left(
0,
a_i+b_{j_s(i)}
\right)
です。
したがって二乗項として
\frac{
r_i^2(s)
}{
a_i+b_{j_s(i)}
}
が出てきます。
ここだけを見ると、同時尤度から導いた式と全く同じです。
実際、
Q(s)
=
\sum_i
\frac{
r_i^2(s)
}{
a_i+b_{j_s(i)}
}
という二乗部分は一致しています。
では、
残差から考える方法と、元の2つの観測の同時尤度から考える方法は全く同じなのか？
というと、ここは少し注意が必要です。
5. 「同じ二乗項」と「同じ尤度」は別
残差 $r_i$ そのものをGaussian変数として正規化まで含めると、
-2\ln L_r(s)
=
\frac{
r_i^2(s)
}{
a_i+b_{j_s(i)}
}
+
\ln
\left[
2\pi
\left(
a_i+b_{j_s(i)}
\right)
\right].
となります。
一方、元の2観測
f_i
\sim
\mathcal{N}(m_i,a_i),
\qquad
g_{j_s(i)}
\sim
\mathcal{N}
\left(
m_i,b_{j_s(i)}
\right)
の同時尤度を書いて $m_i$ を最適化すると、
-2\ln L_{\mathrm{prof},i}
=
\frac{
[f_i-g_{j_s(i)}]^2
}{
a_i+b_{j_s(i)}
}
+
\ln(2\pi a_i)
+
\ln
\left[
2\pi b_{j_s(i)}
\right].
となります。
したがって、
\boxed{
\text{二乗項は同じ}
}
ですが、
\boxed{
\text{正規化項まで含めた尤度は一般には同じではない}
}
ことが分かります。
この違いは、何を「観測データ」として尤度を作ったかに由来します。
6. なぜ違うのか
元の観測は $f_i$ と $g_j$ の2つです。
この2つには、
- 互いの差がどれくらいか
- 2つの共通した強度がどれくらいか
という2種類の情報があります。
残差
r_i=f_i-g_j
だけを見ると、そのうち「差」の情報だけを取り出しています。
一方、元の $f_i$ と $g_j$ の同時尤度では、2つの観測値をそのまま保持した上で、共通強度 $m_i$ を未知量として扱っています。
したがって、
残差だけを確率変数として扱う問題
と
元の2つの観測を同時に扱い、共通強度について最適化する問題
は、二乗項は同じでも、統計的には完全に同じ問題ではありません。
実際、残差の尤度は、共通強度 $m_i$ をある意味で積分して消去する考え方に近く、一方で上の導出は $m_i$ を最適化して消去しています。
この違いが正規化項に現れます。
7. では、比較範囲の外に出たbinはどうする？
ここで自分が一番引っかかったのがこの点でした。
$f$ 側の比較範囲を固定して $g$ を動かすと、位置ずれ $s$ に応じて比較に使われる $g_j$ が変わります。
例えば、
s=0:
\qquad
g_1,g_2,g_3
を使っていたものが、
s=1:
\qquad
g_2,g_3,g_4
になるような場合です。
すると、
$g_1$ や $g_4$ は位置ずれによって尤度から出たり入ったりするのか？
という疑問が出ます。
ここで大事なのは、
比較範囲から外れたことと、観測データそのものが消えたことは違う
という点です。
8. 比較相手がない観測も消えるわけではない
例えば、ある位置ずれでは $g_j$ に対応する $f_i$ がないとします。
この $g_j$ は観測された値なので、元の同時尤度の立場ではデータとして残っています。
ただし、対応する別の観測がないので、
g_j
\sim
\mathcal{N}(u_j,b_j)
という、そのbinだけの未知の真値 $u_j$ を持たせることができます。
$u_j$ は完全に自由なので、
\widehat u_j=g_j
となり、
\frac{
(g_j-\widehat u_j)^2
}{
b_j
}
=0.
したがって、そのbinは位置ずれ $s$ に対して何の制約も与えません。
しかし正規化項
\ln(2\pi b_j)
は残っています。
つまり、
比較範囲の外にある観測は、二乗項には寄与しないが、元の観測データとしては尤度に残っている
と考えられます。
9. 全観測を固定すると正規化は定数になる
全観測データを
\mathcal{D}
=
\{
f_1,\ldots,f_N,
g_1,\ldots,g_M
\}
とします。
位置ずれを変えても、この観測データそのものは変わりません。
したがって元の同時Gaussian尤度の正規化は常に
C
=
\sum_{i=1}^{N}
\ln(2\pi a_i)
+
\sum_{j=1}^{M}
\ln(2\pi b_j).
です。
これは $s$ に依存しません。
その結果、
-2\ln L_{\mathrm{prof}}(s)
=
Q(s)+C
となり、
Q(s)
=
\sum_i
\frac{
[f_i-g_{j_s(i)}]^2
}{
a_i+b_{j_s(i)}
}
だけを最小化すればよいことになります。
この意味で、通常使われているKatsuda型の統計量は、
元の2つの観測プロファイルを同時Gaussian尤度として扱い、未知の真のプロファイル強度を最適化して消去したときの、位置ずれに依存する部分
と解釈できます。
10. では残差から考える方法は間違いなのか？
そういうわけではありません。
固定された1組について
r=f-g
を考えれば、
\mathrm{Var}(r)
=
\sigma_f^2+\sigma_g^2
なので、Katsuda型の二乗項を直感的に理解できます。
したがって、
なぜ分母が $\sigma_f^2+\sigma_g^2$ なのか？
を理解するだけなら、残差から考える方法はとても分かりやすいです。
一方で、
- 位置ずれを変えると使われるbinが変わる
- 比較範囲外の観測をどう扱うか
- 正規化項をどう考えるか
まで含めて尤度として整理したい場合には、最初から元の $f$ と $g$ の同時尤度を書く方が問題設定が明確になります。
つまり、
\boxed{
\text{残差による説明}
\quad\rightarrow\quad
Q(s)\text{ の直感的な理解}
}
に対して、
\boxed{
\text{元の観測の同時尤度}
\quad\rightarrow\quad
Q(s)\text{ の統計的な位置づけ}
}
と考えると、自分には整理しやすく感じました。
11. 一つ重要な注意
比較範囲外の観測は未知の真値によって完全に説明できるので、その二乗項は0になります。
そのため、位置ずれによって対応するbinの数そのものが変わる場合には注意が必要です。
例えば、大きくずらすほど重なる領域が減る設定では、
Q(s)
に含まれる制約そのものが減り、単純に小さい値を取りやすくなる可能性があります。
したがって実際の解析では、探索する位置ずれの範囲全体について共通して比較できる領域を使うなど、
各trial shiftで同じ数・同じ種類の制約を比較する
ことが重要になります。
まとめ
以前の記事では、
Q(s)
=
\sum_i
\frac{
[f_i-g_i(s)]^2
}{
\sigma_{f,i}^2+\sigma_{g,i}^2(s)
}
というKatsuda型の統計量を、「両方に誤差があるのでvarianceを足す」という観点から説明しました。
この説明自体は二乗項を理解する上で自然です。
ただ、あらためて尤度から考えてみると、
f_i
\sim
\mathcal{N}(m_i,\sigma_{f,i}^2),
g_{j_s(i)}
\sim
\mathcal{N}
\left(
m_i,
\sigma_{g,j_s(i)}^2
\right)
という、2つの観測に対する同時Gaussianモデルから出発することで、Katsuda型の式をより対称的に導くことができます。
また今回、自分が一番理解できていなかったのは、
比較範囲から外れた観測は、統計的に消えるわけではない
という点でした。
比較相手がない観測は、自由な未知強度によって完全に説明されるため二乗項への寄与は0になります。一方、Gaussianの正規化項は元の全観測について残るので、全体としては位置ずれに依存しない定数になります。
なお、残差だけをGaussian変数として正規化まで含めた尤度と、元の2つの観測の同時尤度を共通強度について最適化したものは、同じ二乗項 $Q$ を持ちますが、一般には同じ尤度ではありません。
この区別まで含めて考えると、普段使っているprofile-shiftのカイ二乗が「何を仮定した統計量なのか」が、以前より少し明確になりました。