誤差について
撮像データ解析周りの誤差について、定式化の仕方と実例をいくつか解説する。

八木雅文

誤差とは？
真の値からの測定値のズレ
ここで「真の値」とは何かが問題。
測定値から、どのような値を求めたいと考えているのか？
ちなみに測定値や推定値からどの範囲に真の値があるかの指標は「不確かさ」

系統誤差とランダム誤差
系統誤差：
何回やっても同じ量ずれる誤差
ランダム誤差：
毎回再現性のない誤差
例えば機器の目盛りが間違っていて毎回同じ量ずれる⇒系統誤差
例えば機器が不安定で結果に再現性がない⇒ランダム誤差

真の値？
誤差は真の値からの測定値のズレ変化しないものは、真の値を言うのは簡単。
例えば真空中の光速度は厳密に
毎秒 299792458 m
一方、光の強度の場合、光子はミクロな確率過程を経て天体から出てくるため、ある１秒間に天体から出てくる光の量は次の１秒と同じではない。
つまり観測値は本質的に一定ではない。

光量の誤差とは？
その上で、天文の測光の場合、求めたい「真の値」は、観測される光量の「確率分布の代表値」
確率分布が時間変化しない場合には
これは「長時間平均」であると考えてよい。
光量の誤差は、測定あるいは推定した値がどれだけ精度良くその確率分布の代表値を求められたかという指標。

たとえ話
東京都の1日の交通事故数を調べたいと考え、ある日が150件、次の日は140件という情報を得たとする。
この生データ値を可能な限り正しく求めるのも重要だが、最終的に欲しいのは「典型的な発生件数」なので、150とか140という数字自体の精度だけをひたすら追求するのはあまり意味がない。
欲しいのは、例えば 「平均 142.0 件/日」という情報。

たとえ話の蛇足
直前のスライドは2011年版の再利用なので例が交通事故でコロナではないのだった、
とかいう背景はさておき、
コロナ陽性者数の場合は時間変化する。
このような時間変化する量の場合も、欲しいのはその時の厳密な発生数ではなく、その時の確率分布の典型値。
天文の場合、いくつ光子が来たかではなく来るはずだったか、を知りたい。

目指すところ
光量の推定の場合、光子数のポアソン誤差は原理的に減らすことができないので、目指すのは、ポアソン誤差以外の誤差を
ポアソン誤差に比べて十分無視できるくらいまで小さくすること。
なお、HSCの観測では多くの場合、
目標天体以外の天体、主に夜空からの光子のポアソン誤差が一番大きい。
空が暗い時に観測したい理由の１つはこれ。

確率分布観測量Ｘ：観測ごとに変わる量があったとして、ある試行でx0である確率を p(x0)と書ける場合、
関数p(x)が確率関数。
Ｘがとる値xが連続量の時、
x0-dx/2<x<x0+dx/2 に入る確率がp(x0)dx になるようなp(x)が、確率密度関数。
これらを本資料中では(確率)分布と呼ぶことにする。

分布と代表値分布の代表値として良く使われるのは、
平均、中央値（メディアン）、最頻値（モード）など。
今回は平均を代表値と考える。
頻度

値
平均
最頻値
中央値

誤差の評価一回の測定での誤差は、真の値がわからないと求まらないが、測定値の「分布」が仮定されていたり何らかの方法でわかっていたりすれば「誤差の典型値」を計算できる。
例えば良く使われるのが標準偏差（誤差の二乗平均）σなどと書かれる。
※σは正規分布の標準偏差に相当する別の量の場合もある。

誤差の推定１）そもそも誤差はどうなるはずかから考え、誤差の伝播を使って観測量の誤差を推定する(測光の限界等級)２）観測量の統計を考えた上で、解析での誤差の広がりを推定する(画像などの並行移動)を例として挙げる。

二項分布
ある確率pで事象が起きる試行をＮ回行ったとき、Ｘ回事象が起きる確率分布を二項分布と呼ぶ。書き方には色々あるが、ここではB(N,p) と書くことにする。
x
N −x
B ( N , p ) ( x)= N C x p (1− p)
N!
NC x =x ! ( N − x) !

例えば N=5, p=1/6 の時、(６面サイコロを５個振った時に１の目の出る数の確率分布）はB(5,1/6) と書ける。B(5,1/6)(0)～0.402
二項分布の平均は Np, 分散は Np(1-p)

ポアソン分布
二項分布で Np を一定のまま p→0 Ｎ→∞に飛ばした極限がポアソン分布。
k=Np と置いて、p(k)は以下になる。
平均 E(X)=k
分散 V(X)=k

X

−k

k e
p (k )( x )=x!

二項分布では平均は Np、
分散は Np(1-p)だったが、ポアソン分布は p→0 の極限なので、1-p→1でk=Np=Np となる。
標準偏差は√ｋポアソン分布する変数の和
ポアソン分布に従うＸ、Ｙという２つの変数があったとき、和 Ｘ＋Ｙ もポアソン分布に従う。証明略。
一方、二項分布の和は、ｐが等しければ二項分布になるが、pが違うと二項分布には必ずしもならない。
また、差、Ｘ－Ｙはポアソン分布には従わないので要注意。

ポアソン分布する変数の定数倍
ポアソン分布に従う変数の定数倍はポアソン分布に従わない。
ここも要注意。
電荷数はポアソン分布に従うが、アンプを通して出てきた値はゲインで割り算されているので、ポアソン分布には従っていない。
変数に定数を足してもポアソン分布に従わなくなる。

ポアソンを確率で間引く
例えば単位時間当たりの個数がポアソン分布に従う光子が、透過率 f の系を透過して来た結果もポアソン分布に従う
二項分布にまで戻って考えればB(N,p) の確率部分pが p で発生した後 f の確率で生き残るので確率が fpになったのと同じ。
ポアソンでも期待値kや分散kがfkになる。

ポアソンを確率で間引く
解析的にも確認しておく。間引かれた結果zを観測し、間引かれたのがｙだった場合、( y + z) −k
その確率は ( y + z) !
k
e
z

y! z!

これを足していけば( y +z ) −k∞

∑ y= 0kp (1− p)z −ky( y + z )!

(kp) e(k (1− p))p (1− p) =
∑
y=0
z!
y!
z −k
z −kp
( kp) ek (1− p) (kp ) e
.=
e
=
z!
z!

e
y! z!

z

y

∞

y

と、平均 kp のポアソン分布にちゃんとなる。

正規分布
別名ガウス分布。
平均ｍ、標準偏差σの場合、

 x−m 
−
2 2

2

1
N  m ,  x=e2
2  
2

ポアソン分布は「kが大きければ」
正規分布 N(k,k) で良く近似できる。
多くの誤差は正規分布だと仮定されることが多い。

kが大きければ？
どれだけkが大きければ、正規分布は
ポアソン分布の良い近似になるのか？
これ、何をもって「良い近似」とするか、
色々考え方ありますが、k=10と100をグラフにすると、こんな感じ。

誤差の計算
確率変数、Ｘ，Ｙと関数f(X,Y）があった場合、f(X,Y)の誤差は概ね以下の方針で求められる。
１）確率変数Ｘと誤差δＸを求める。
Ｙも同様にδＹを求める。
２）期待値E(f(X,Y))を求める
３）分散V(f(X,Y))=E(（f(X,Y)-E(f(X,Y))^2)
つまり観測値から期待値を引いた二乗平均を求める。
この平方根が誤差の典型値。

誤差の伝播の例
Ｘ，Ｙの平均を  X , Y と書き、標準偏差を  X ,  Y とする。
この時、 E ( X +Y )=μ X +μY
2

V ( X + Y )=E ((( X + Y )−(μ X + μY )) )2
=E  X − X Y − Y  
2
=E  X − X  2 E  X − X Y − Y 
2
E Y − Y  22=σ X +σ Y +2 E (( X − μ X )(Y − μ Y ))なお、一般的にE(aX+bY)=aE(X)+bE(Y)共分散
ここで出てきた E (( X − μ X )(Y − μ Y ))を共分散と呼んで Cov ( X , Y ) と書く。
共分散が０の時、Ｘ、Ｙは独立という。
2
2
この時は、V ( X +Y )=σ x +σ y
誤差の伝播で変数が独立かどうかは極めて重要な点なのだが、実際のデータでは独立かどうか判断するのは、必ずしも簡単ではない。
22x

追記：あとこれも有用。 V (aX )=a σVが定義できる分布なら常に成立する。

バイアスの統計
例えばバイアスを補正する場合、データ１枚のバイアスが N(b0,r^2),
バイアスパターンはn枚合成してN(b0,r^2/n) となっている場合、独立ならば差し引いた後の統計は N(0,r^2(1+1/n))
読み出しノイズ相当の量は sqrt(1+1/n) 倍だけ増えている。n=10 の場合、5%増。
読み出しノイズを5%増やしてでもバイアスを引いたほうがいいのかは考える必要がある。
（なお通常はポアソンノイズが圧倒するので読み出しノイズが5%増えても気にならない）

追記：平均と分散
E (aX +aY )=a μ X + bμY
2
2
V ( X +Y )=σ X + σ Y + 2Cov ( X ,Y )2 2
V (aX )=a σ xなどは、平均や分散が定義できれば、X,Yが正規分布でなくても成立する
強い式なので使いでがある。
「定義できれば」とか書いているのは、
コーシー分布のように平均や分散が定義できない分布もあるため。

限界等級
撮像データの測光の限界等級とは、天体の測光値（Signal)が、測光誤差(Noise)の何倍かという比（S/N） がある一定値、例えば、10以上の天体が何等か、という値。
例えば28等の天体が S/N=10 の時、「S/N=10 での限界等級が28等」
と言う。

限界等級
この場合の限界等級は「測光精度」の意味での限界等級であって、検出の限界の意味での限界等級ではない。
測光精度の意味での限界等級が同じであっても、天体検出の割合が同じであるとは全く限らない。
全然別の話である。

検出限界
同じ等級でも星のような点光源では、ほとんどの光が小さな範囲、例えば
半径5pixくらいに収まるのに対し、銀河ではもっと広い範囲に光が淡く広がる。
このため、銀河の方が検出が難しく、特に大きく広がったものの検出はどんどん難しくなる。
（測光のS/Nも悪くなる）

検出限界
検出限界は私の知る限りでは、解析的に求めるのは極めて困難。
（もし解析モデル知ってる方
おられましたら御教示ください）
多分、ある確率で画素が閾値を超えるような画像で、閾値を超えたピクセルがN個つながる
確率は？、という問題になる、はず。
実用上は、疑似画像を埋め込んで検出確率を求めることが多い、と思う。

等級の系統誤差とS/N
等級の場合、まず単位時間あたりの光子数がポアソン分布し、フラックスは正規分布と仮定できる。このときフラックスを F, 等級をmと置くと、等級の期待値はE (m)= E (−2.5 log( F )+ const)
2.5
E (m)=−
E (ln ( f 0 +δ F ))+ constln (10)δFE (ln ( f 0 +δ F ))= E (ln ( f 0 )+ln ((1+
)))
f
0
2
3
1 δF
δ
F
δ
F
δ
F
ln 1+
=
−
+O
f0
f0 2 f0
f0()( ) ( )

等級の系統誤差とS/N
Fは正規分布だと思っているのでδFの平均は02
2
E ( ( δ F ) ) =σ((δFと置けば、 E ln 1+ f
0

))

(

1 σ 2
1
1
=−
=−
2 f0
2 (S / N )( )したがって、 m0 =−2.5 log ( f 0 )+constとしたとき、
22.5
1
0.54
E (m)=m0 +∼m 0 +
2
2 ln (10 ) ( S / N )(S / N )()

S/Nが7くらいより小さいと、0.01等以上の「系統誤差」が等級に乗る、と計算できる。

)

2

等級の系統誤差とS/N
・・・書いていて本当にそうなるのか不安になってきたので一応確認してみた。
S/N=5 の場合を想定して、平均５標準偏差１、N(5,1^2)の乱数を10000回振る。
その結果を -2.5logして求めた等級の平均をとって、-2.5log(5) と比較してみると、
確かに 0.022 等暗い側にズレている。
0.54/5^2=0.022 なので、計算は合っている。
なお、このトイシミュレーションでは等級のメディアンは 0.001 等しかずれなかったが理論的にもそうなるのかは未確認。

等級のランダム誤差
一方、等級のランダム誤差はE((m-m0)^2) で評価できる。

(( ) ) なので、

2.5 δ FδF 2
m∼m 0−+Oln (10 ) f 0
f0

(

2.5
V (m)∼ln (10))

2

1( S / N )2(

σ (m)∼

)

2.5
1
1.086
∼
ln (10) S / N S / N

例えば 0.1 等以下の精度の測光をしたいのであれば、S/N>10.86 が必要。

等級の系統誤差とS/N
S/N=5 のときNに正規分布乱数を100000回振って書いた等級誤差のヒストグラム

足し合わせの方法
標準偏差で誤差を評価した場合、S/N を最も高くする足し合わせ方は重みづけ平均。
しかし、重みづけ平均は外れ値に弱い。
そこで、複数の露出の足し合わせにはsigma-clipped mean や、メディアンが用いられる。

メディアン
例えば、平均１００、分散20の観測値のうち１つに非常に大きなノイズが乗った場合、
113, 92, 66, 10000, 104, 68, 90
これを単純平均を取ると ～1500。
一方、median は 92,
MAD でσを推定した上での、
3σ-clipped mean は、89
単純平均よりは遥かに良い推定

メディアンの誤差
しかし、メディアンは単純平均よりも誤差が大きくなる。
分散１の正規分布の平均とメディアンを比べると、Ｎが十分大きいところでは、
標本標準偏差は約25%悪化する。
1s (mean)=√N
1.25
s median=Nメディアンの誤差
Ｎ s(median) s(1-clipped)
3 1.161.22
4 1.091.15
5 1.201.12
6 1.141.10
7 1.211.08
単純平均の標準偏差(s(mean))との比。
この値が大きいほど誤差が大きい。
１つだけclip される場合とメディアンではN>4 では clipped の方が誤差が小さい。
なお偶数 median で誤差が小さいのは
中央２つの平均を取っているから。

系統誤差の影響
１σ程度の系統誤差、例えば画像の足し合わせの際のバックグラウンドの引き残りやゴーストなどには、メディアンは引きずられにくいが、
clipped mean は clip がうまくできないと引きずられる。
一方、外れ値が常に大きい側に偏る
場合などは、メディアンは系統的にそちら側にずれる。

メディアンのずれ
具体的には例えば N=2m+1 個のデータのメディアンで、1個大きい側に外れ値が混じる場合、本来は大きい方から m+1 番目をとるはずが、m 番目をとることになる。
このことでmedian xは、元のデータが外れ値以外は N(μ,σ^2)に従う場合、
p(x)=(2m)!/((m-1)!(m+1)!)(1-erf(x))^(m-1)gauss(x)erf(x)^(m+1)に従うことになる。

メディアンのずれ
定量的には、元のデータが外れ値以外 N(μ,σ^2)に従うとして、
サンプル数Nが十分大きく、大きい側の外れ値の混じる割合が h (h<0.1程度)の場合、経験的には～1.25hσほどプラス側にズレるようである。
Nが小さい場合、1つ大きな外れ値が混じった場合、N=3,5,7 では
0.5σ,0.3σ,0.2σほど、平均値がプラス側にずれる。

メディアンのずれ
外れ値が１つ入った場合の分布は
こんな感じになる。N=5 で、１つ外れ値が混じると、0.3σほどずれる。
この分布のｘの一次のモーメントが外れ値入りの場合のメディアンの期待値（平均値）

メディアン
したがって、常に外れ値がプラスにしか入らないような場合には、
メディアンは系統的にずれる。
これを考えると、「もしうまくクリップができるのであれば」
クリップ平均の方が結果が偏らない。

追記：やってみた
講義後、宇宙線が入ったらプラスに偏るのか？
という質問をいただいたので「そうです」と答えつつ念のため実際に試してみた。
bias ５枚からあるチップあるチャネルの適当な範囲を切り出して、そのうち１枚のある領域に+10000
して、残り４枚とmedian stack した結果３種。

四角はわかる。幅１の縦線はわからない（と思う）
幅３の斜め線はかろうじてわかるか。
と、目で見てわからなくても、このように偏る。

画像の並行移動
回転等でも類似の問題があるのだが、ここでは並行移動を考える。画像・データを非整数画素だけ並行移動させたい。
例えば位置を合わせて足し合わせる場合など。
方法として、線形補間のリサンプリングを考える。
例えばｘ軸方向に -a (0<a<1) シフトさせるときはv ' ( x ' )=(1−a ) v ( x)+a v ( x +1)

x x+a x+1

この結果、１画素あたりの分散は ((1− a)2 + a 2 ) σ 222=(1−2 a (1−a )) σ <σと小さくなる。

実際にシフトさせてみた
左がシフト前、右はx+0.5, y+0.5 シフト

見た目でもざらつきが減っているのだが、実際測ってみると、画素あたりのノイズは確かに減っている。左67ADU, 右43ADU
すると並行移動でS/Nが上げられる！？

もちろんそんなことはない
シフト後は画素間が独立ではなくなる。
このため、ある範囲の画素を足した合計の分散は、各画素での誤差の二乗（分散）の和に加えて、更に共分散を足す必要がある。
具体的に[0,N-1] の範囲の総和の誤差を計算してみる。
元の誤差は独立なN点の値の和なので分散はV (v (0)+v (1)+... v ( N −1))= N σ

平行移動後

2

v ' (0)=(1−a ) v (0)+a v (1)v ' (1)=(1−a ) v (1)+a v (2 )...

誤差は減っていない
V (v ' (0)+v ' (1)+... v ' ( N −1))
=V ((1−a ) v (0)+v (1)+... v ( N −1)+ a v ( N ))2222=((1−a ) +( N −1)+a ) σ =( N −2 a (1−a ))σ

Nが大きければほとんど元と変わらないまま。
このような１画素あたりの分散が見かけ上小さくなることは画像処理の結果よく起きる。
このため、実際のデータの測光誤差の推定は画素ごとではなくある程度の範囲のカウントの和のバラツキで通常は推定する。

観測モデル
モデルとして以下を考える。
( OBJECT +SKY )× flat × t obs +DARK × t dark
DATA =+ BIAS
gain

OBJECT,SKY,DARK,BIAS は確率変数で、観測量の DATA も確率変数。
flat, tobs, tdark, gain は定数。
モデルとしては１回の観測では分布に従った
確率変数 OBJECT,SKY,DARK,BIAS から
DATA が得られると考える。

仮定
恐らくは、画素１つ当たりの光電荷の数はポアソン分布に従い、そこを読み出した
カウントは正規分布（ガウス分布）に従うだろう、と「仮定」する。
独立性に関しては、あるピクセルでの光電荷の数は隣のピクセルの光電荷の数とは、独立だろうと仮定する。
※ただし、実際それぞれの装置で、画素間の独立性がどこまで成り立っているかは、
私は未確認です。

S/N のモデル計算
露出時間をt、システムの透過率をf、天体の光度をLとすると、検出器で検出される信号は( OBJECT +SKY )× flat × t obs +DARK × t dark
DATA =+ BIAS
gainのうち OBJECT × flat ×t obs の部分。
単位は光電荷数、としておく。
単位時間当たりの光子数（相当の量）に時間と透過率を掛け算。

追記：flatの誤差講義で、flat の誤差について質問いただき、「実際にはflatの誤差はすごく効く」と答えたのですが、その辺りを追加解説。前のページの式はあくまで我々は実際には知らないけど、理想的にはある定数としての flat(透過率),gain があってそれは試行を繰り返しても変わらないので確率変数ではなくモデル中では誤差もないというモデルです。強いて書くなら真の値 flat0, gain0とでも書くべき値でした。
実際の解析では観測データから推定した
flat=flat0+Δflat を使うことになるのですが、この
flat の誤差Δflat がDATAからOBJECTを求める際にすごく効くはずです。ただし、全く同じ条件で標準天体で較正できれば、この誤差はキャンセルして消えます。

追記：flatの誤差
横道にそれますが flat の誤差Δflat は、
ドームフラットの各ピクセルのポアソン誤差、
フラット板や光源の非一様性、観測時の姿勢と違うことによる周辺減光パターンの違い、
ほとんど効かないけど読み出し雑音やダークの誤差に加えて、全体をcoaddして足し合わせる前の各データの規格化の誤差などなど、考えるべき要素はあれこれとあるので結構大変です。
しかも多分1%未満のレベルまで気にして調査する必要があるのでそれも大変です。

S/N のモデル計算
OBJECT × flat ×t obsのOBJECT はポアソン分布に従う変数
式は定数を掛けているように見えるが、これは変数の定数倍ではない。
flat は、確率での間引き、
tobs は最初の１秒の光量のtobs倍ではなく、整数で考えれば独立な
tobs 回の光量の和。
OBJECT の期待値を object と書いてS = object × flat × t obs

S/N のモデル計算
次に、S/N の N の方。誤差。
測光での大きな誤差源の一つは機器のノイズ（読み出し雑音）。読み出し雑音は
露出時間によらず一定であると、とりあえず考えてよい。
１画素あたりのこの誤差を ｒ として、
E ( N 1 )=0 ,V ( N 1)= r

2

天体が A 画素に入っている場合、読み出し雑音は画素ごと独立だと仮定してE ( N ' 1)=0 ,V ( N ' 1 )= A rと、分散はA倍。

2

S/N のモデル計算
もう一つの大きな誤差要因は測定される電荷の確率的な誤差。
（ポアソン雑音）
目標天体の画素数をＡとすると、ポアソンの期待値は、OBJECT,DARKの期待値を object, dark と書いて、(object +sky × A )× flat ×t obs + dark × A × t darkとなる。分散（誤差の二乗）は、V ( N 2 )=( object + sky × A )× flat ×t obs + dark × A × t dark

測光のS/N
読み出し雑音とポアソン雑音は全て独立と考えられるので、雑音の分散は合わせて2

V ( N )= Ar +( object + sky × A )× flat ×t obs +dark × A ×t dark

誤差NはN =√ Ar 2 +(object + sky × A )×flat × t obs + dark × A ×t dark
object ×flat × t obs
S =N √ Ar 2 +(object + sky × A )×flat × t obs + dark × A ×t dark

測光のS/N
object ×flat × t obs
S =N √ Ar 2 +(object + sky × A )×flat × t obs + dark × A ×t darkｔが小さい場合、分母の平方根の中は Arが残って、 S ∼ object ×flat × t obs
Nr √A

ｔが大きい場合、 t obs ∼ t dark として、
S ∼object ×flat
× √ t obs
N √ ( object +sky × A )× flat + dark × A

2

つまり

S
露出時間tが小さい場合は、 N ∝t
S
大きい場合は、 N ∝ √ t

露出時間が長い場合、S/N は
露出時間の平方根に比例。
S/N を倍にしたい場合は
露出時間は４倍必要、精度を１桁上げるには１００倍必要。

まとめ
このようにモデルを立てて、順に計算すれば誤差がどうなるか推定できるケースがそれなりにある。
しかし、一方で解析的な推定が難しい問題もまた多く存在する。
そのような場合、例えば検出限界のように
疑似データを加えて実際に検出を試してその結果から推定するなど、シミュレーションから誤差を求める場合も多い。