考える波 第 59 回 / 第 XIX 部・道具を持って 物理へ戻る

系統誤差を「見積もる」から「測る」へ

系統誤差を測る ──
そしていちばん効いたのは、窓ではなかった 第 58 回の「窓を変えると値が動く」を逆手に取り、動かない区間を探す手続きを作りました。
★★★ 途中で分かったこと ── あの「行きすぎ」は雑音ゼロでも出ます。主役は記録の分解能より低い周波数、つまり記録の中では傾きにしか見えない成分でした。
★★ だから直線を 1 本引くだけで消えます(\(\beta=5\) で \(5.512\to4.992\))── 第 58 回で「窓の階数を上げろ」と書きましたが、もっと安い手がありました。

必要な道具:第 51, 55, 57, 58 回計算:kenshou/keito.js

第 58 回で符号の決まった系統誤差が出ました。★ 今回はそれを測定の道具に変えます ── 動かし方から、真の値と誤差の両方を取り出す。

01\(\hat\beta\) を \(\alpha\) について走らせる

窓の階数 \(\alpha=1\)(矩形)から 8 まで。32 realization の平均です。

真の \(\beta\)\(a{=}1\)2345678
10.9900.9890.9910.9880.9880.9870.9880.990
32.0333.1602.9922.9902.9892.9882.9892.990
51.9844.0205.5125.0614.9904.9904.9904.991
61.9814.0216.0007.6235.9975.9925.9905.992
71.9814.0216.0048.0017.6337.1826.9926.994

01節の中心

★ 小さい \(\alpha\) では漏れで頭打ち(第 58 回の天井 \(2\alpha\))。
★★ 天井のすぐ手前では行きすぎる ── \(\beta=6\) で \(\alpha=4\) が 7.623、\(\beta=7\) で \(\alpha=4\) が 8.001。いずれも天井 \(2\alpha=8\) に向かって突き上げています。
大きい \(\alpha\) では落ち着き、どの \(\beta\) でも真値より 0.01 ほど低い一定値に並びます ── ここが「平らな区間」。
★★★ 手続きの骨はこれです平らな区間を自動で見つけ、その幅を系統誤差として報告する。

02★★★ 雑音を消して、「行きすぎ」が本物か調べる

期待される周期図 \(E[P(k)]=\int S(f)|W(k-f)|^2df\) を直接たたみこんで、統計誤差ゼロの偏りだけを出します。

真の \(\beta\)\(a{=}1\)2345678
10.0050.0000.0000.0000.0000.0000.0000.000
3−1.009+0.1250.0010.0010.0010.0010.0010.002
5−3.012−0.988+0.625+0.0680.0030.0030.0040.004
7−5.012−2.988−1.000+1.002+0.918+0.2520.0070.008

★★★ 行きすぎは雑音ゼロでも出ました(\(\beta=5\)、\(\alpha=3\) で \(+0.625\))── 第 58 回で測った \(+0.61\) とほぼ一致。統計の効果ではなく、期待されるスペクトルそのものの形です。

★★ 予想が外れました大きい \(\alpha\) で主ローブの平滑化による偏りが育つと思っていましたが、\(+0.002\sim+0.008\) しかありません。大きい \(\alpha\) 側を止めているのは偏りではなく、ばらつきでした(04 節)。

では、どこから漏れているのか。積分の下端を動かしてみます(\(\beta=5\)、\(\alpha=3\))。

積分の下端 \(f\)[bin]0.06250.250.51.0
推定 \(\beta\)5.6255.0625.0065.002

02節の中心 ── 漏れの主役が分かった

★★★ 記録の分解能 1 bin より下(0.0625〜1 bin)を入れたときだけ、行きすぎが出ます。
★★ 漏れの主役は「記録の中では傾きにしか見えない成分」でした ── 周期が記録より長い波。第 5 回・第 56 回で「観測時間より長い周期は定数と区別がつかない」と書いた、まさにその成分です。
区別がつかないだけでなく、測定を汚していた

03★★★ ならば、傾きを引けばよい

記録の中で傾きにしか見えないなら、当てはめて引いてしまえば消えるはずです。ハン窓(\(\alpha=3\))で、\(q\) 次多項式を引いてから測ります。

真の \(\beta\)引かない\(q{=}0\)(平均)\(q{=}1\)(直線)\(q{=}2\)\(q{=}3\)
32.9922.9922.9912.9912.991
44.0004.0003.9923.9923.991
55.5125.5124.9924.9934.992
66.0006.0006.0355.9945.992
76.0046.0046.5046.0396.574

03節の中心 ── 本回いちばんの実用的な収穫

直線を 1 本引くだけで、行きすぎが消える。
\(\beta=5\):\(5.512\to4.992\)

★★★ 第 58 回では「窓の階数を上げろ」と書きましたが、もっと安い手がありました。階数を上げるとばらつきが増えます(04 節)が、トレンドを引くのはほぼ無料です。
平均を引くだけ(\(q=0\))では効きません ── 効くのは 1 次から。
★★ ただし天井は動きません ── \(\beta=7\) を \(\alpha=3\)(天井 6)で測ると、どの \(q\) でも 6.04〜6.57 のまま。トレンド除去は「行きすぎ」を消しますが、「貼りつき」は消しません。
引きすぎも害になります(\(\beta=7\) で \(q=3\) にすると 6.57 と悪化)。手引きに足すべきは「1〜2 次を引き、そのうえで \(2\alpha\ge\beta+2\) になるよう窓を選ぶ」。

04ばらつきは \(\alpha\) とともに増える

窓をかけると独立な点が減ります。等価雑音帯域 \(\mathrm{ENBW}=N\sum w^2/(\sum w)^2\)。

\(\alpha\)ENBW[bin]実測 sd(\(\beta\))第 51 回の予言
11.0000.11700.0865(偏りあり)
21.2340.18150.0961(偏りあり)
31.5000.08390.10590.792
51.9440.10110.12060.838
82.4730.11930.13600.877

★ 予言は第 51 回の \(\mathrm{sd}=1.2825/\sqrt{\text{有効点数}}\)、有効点数 \(=S_{xx}/\mathrm{ENBW}\)。\(\alpha=1,2\) は \(\beta=3\) では偏っているので、この sd は意味を持ちません。

★★ 偏りの無い \(\alpha\ge3\) では比が 0.79〜0.88 ── 予言は 15〜20 % 大きめでした。周期図の点が窓のせいで相関しており、ENBW だけではその相関を過大に見積もるためです。使うなら上限として。

05最良の \(\alpha\) ── ただし乗ってはいけない最小値がある

真の \(\beta\)\(a{=}1\)2345678
31.0110.2420.0850.0930.1000.1060.1120.116
42.0190.0260.0850.0930.0990.1050.1110.116
64.0191.9780.0081.6860.0970.1040.1100.115
75.0192.9790.9981.0010.8850.2540.1100.115

05節の中心(表は RMS 誤差)

★★★ \(\beta=4\) の \(\alpha=2\)(0.026)や \(\beta=6\) の \(\alpha=3\)(0.008)は、偏りがたまたま打ち消し合っただけ。\(\beta\) を知らなければ狙って当てられません ── この最小値に乗ってはいけない。
隣を見れば分かります:\(\beta=6\) で \(\alpha=3\) は 0.008 ですが、\(\alpha=4\) は 1.686。これほど不安定な最小値は、使える最小値ではありません。
★★ 安全な目安は \(2\alpha\ge\beta+2\)(\(\beta=4\) なら \(\alpha\ge3\)、\(\beta=6\) なら \(\alpha\ge4\))。\(\beta\) を知らないと選べないので、01 節の走査が要ります。

06★★ 手続きにする ── そして当たり具合を測る

手順

① \(\alpha=1\ldots8\) で \(\hat\beta\) を測る ② 幅が \(\mathrm{tol}=0.15\) 以内で続く いちばん長い区間を探す ③ 値はその平均、系統誤差はその半値幅 ④ 統計誤差は第 51 回の式(ENBW 込み) ⑤ 二乗和で足して報告する

真の \(\beta\)平均の推定合計誤差\(|\)誤差\(|\) 中央68 % 点95 % 点実測 sd被覆率
11.0150.1290.0570.0780.1780.31787.0 %
33.0210.1290.0600.0920.1870.36684.0 %
55.0140.1320.0620.0930.2010.22484.5 %
77.0090.1360.0680.1050.2160.12179.5 %

06節の中心 ── 誤差の分布が正規分布ではない

名目 68 % に対して実測の被覆率は 79〜87 % ── 誤差を 1.4 倍 ほど大きめに出しています(\(|\)誤差\(|\) の 68 % 点が 0.08〜0.11 なのに、報告する合計誤差が 0.13)。
★★★ ところが実測 sd は合計誤差よりずっと大きい(\(\beta=1\) で 0.32 対 0.13)。矛盾ではありません ── 誤差の分布が正規分布ではないから。
★★ 芯は細く、まれに大きく外します(95 % 点は 68 % 点の 2 倍 ほど、sd は 0.12〜0.41)。その稀な外れが sd を膨らませています。
★★★ だからこの手続きの誤差は sd で書いてはいけません。「68 % がこの幅に入る」と被覆率で書くのが正しい ── 第 51 回の \(\ln(\chi^2_2/2)\) が非対称だったのと、同じ理由です。

07手続きが破れるところ

当てはめ帯域真の \(\beta\)推定合計誤差被覆率
\(k=32\ldots512\)33.0210.12984.0 %
\(k=32\ldots512\)66.0060.13383.0 %
\(k=64\ldots256\)33.0270.35275.5 %
\(k=64\ldots256\)66.0150.36778.5 %
\(k=16\ldots1024\)33.0190.07486.5 %
\(k=16\ldots1024\)66.0080.07585.0 %

帯域を変えると被覆率も動きます(75.5〜86.5 %)── 推定値そのものはほとんど動かない(3.02 前後)のに。

★★ つまりこの手続きも「仮定つき」(第 55 回)仮定は三つ ── ① 当てはめ帯域 ② tol ③ 窓の族。手続きを作ったからといって、仮定が消えたわけではありません。

練習問題

  1. 02 節で「漏れの主役は記録より長い周期」と分かった。第 56 回とどうつながるか。
    答えを見る
    第 56 回は「見えない」と書きましたが、今回は「見えないのに効いている」が分かりました。
    ★ 第 56 回:「観測時間より長い \(\tau\) は折れ曲がりが見えない」。★★ 本回:それでも裾は窓を通って高周波に漏れ、傾きを \(+0.6\) も汚す。
    ★★★ 第 57 回 07 節の練習問題で「\(n\) が小さいほど端の外の情報が遠くまで効く」と書いたのは、これのことでした ── あのときは「原理的には \(\tau_{\max}\) を推定できるかもしれない」と書きましたが、実際に起きるのはまず汚染です。
    順序としては:まず引いて(03 節)、それから測る。推定に使うのはその後の話。
  2. 03 節で「平均を引くだけでは効かない」と出た。なぜ 1 次から効くのか。
    答えを見る
    窓が 0 次の漏れをすでに消しているからです。ハン窓(\(\alpha=3\))は端で値も傾きも 0 なので、一定成分(\(f=0\))の漏れは \(f^{-3}\) で十分に落ちています。
    残るのは \(f\) が 0 でない低周波 ── 記録の中では「直線」に見える成分。★★ 直線を引くと、その成分の大半が消えます。
    ★★★ 一般化すると\(q\) 次を引くのは「記録の中で \(q\) 次多項式に見える低周波」を消すことで、それは実効的に窓の階数を上げるのと同じ働きをします。違いは代金窓の階数を上げるとばらつきが増えますが(04 節)、トレンド除去は自由度を \(q+1\) 個 減らすだけ ── 4096 点なら ほぼ無料です。
  3. 06 節で「sd で書いてはいけない」と出た。ではシリーズの過去の数字は大丈夫か。
    答えを見る
    第 51 回の \(\mathrm{sd}(\beta)=1.2825/\sqrt{M}\) は「そのまま sd」なので問題ありません ── あれは \(\ln(\chi^2_2/2)\) の標準偏差そのもので、手続きの出力ではないからです。
    危ないのは第 55 回で「模型を仮定した値」に分類した 5 件 ── あれらは手続きの出力に近く、分布が正規とは限りません。
    ★★ 本回の教訓を第 55 回の洗い直しに足すなら「模型を仮定した値」には、誤差の形(正規か、芯が細く裾が重いか)も書く
    ★★★ そして本回のように被覆率を実測するのがいちばん確実です ── 誤差の「大きさ」より先に、誤差の「形」を測る。

まとめ 系統誤差は測れる。ただし、いちばん効いたのは窓ではなかった

第 58 回の「窓を変えると値が動く」を逆手に取り、動かない区間を探す手続きを作りました。

\(\hat\beta(\alpha)\) には平らな区間ができます ── 小さい \(\alpha\) では天井 \(2\alpha\) に貼りつき、その手前で行きすぎ、大きい \(\alpha\) で落ち着く。

★★★ そして途中で分かったこと ──
あの「行きすぎ」は雑音ゼロでも出る。
主役は記録の中では傾きにしか見えない成分だった。

★★ 積分の下端を動かすと決着しました ── 分解能 1 bin より下(0.0625〜1 bin)を入れたときだけ、\(+0.625\) の行きすぎが出る。第 5 回・第 56 回の「観測時間より長い周期は定数と区別がつかない」は、「区別がつかないだけでなく、測定を汚していた」でもありました。

★★★ だから直線を 1 本引くだけで消えます(\(\beta=5\) で \(5.512\to4.992\)、\(\beta=6\) は \(q=2\) で 5.994)。第 58 回で「窓の階数を上げろ」と書きましたが、もっと安い手がありました ── 階数を上げるとばらつきが増えるのに対し、トレンド除去は 4096 点なら ほぼ無料です。ただし天井は動きません(\(\beta=7\) を \(\alpha=3\) で測ると どの \(q\) でも 6.04〜6.57)。

★★ 予想が一つ外れました大きい \(\alpha\) で平滑化の偏りが育つと思っていましたが \(+0.002\sim+0.008\) で無視でき、大きい \(\alpha\) 側を止めているのはばらつきでした。そして「RMS が最小の \(\alpha\)」は偏りの偶然の打ち消しで決まることがあり(\(\beta=6\) の \(\alpha=3\) で 0.008、隣の \(\alpha=4\) は 1.686)、乗ってはいけない最小値です。安全な目安は \(2\alpha\ge\beta+2\)。

★★★ 手続きの当たり具合も測りました ── 名目 68 % に対し被覆率 79〜87 %。ところが実測 sd は合計誤差の 1〜3 倍。矛盾ではなく、誤差の分布が正規でないからです ── 芯は細く、まれに大きく外す。だからこの手続きの誤差は sd ではなく被覆率で書くべきで、それは第 51 回の \(\ln(\chi^2_2/2)\) が非対称だったのと同じ理由です。

そして手続き自体も「仮定つき」でした ── 当てはめ帯域・tol・窓の族 の三つ。帯域を変えると被覆率が 75.5〜86.5 % と動きます。手続きを作っても、仮定は消えません。

この文書は「考える波」シリーズ第 59 回です。計算は kenshou/keito.js で行っています(FFT は自前、依存ライブラリ 0)。本回の数値について:\(\hat\beta(\alpha)\) の走査、雑音ゼロの期待周期図からの偏り、積分下端の依存性、トレンド除去の効果、ENBW とばらつき、RMS 誤差、手続きの被覆率と誤差分位点は、すべて本稿で計算したものです。周期図のスペクトル漏れ、等価雑音帯域、周期図を取る前にトレンドを除くべきこと ── これらはいずれもよく知られた実務で本稿の発見ではありません ── 本稿がしたのは、それらを第 58 回の \(\alpha\) という物差しで整理し、行きすぎの出どころを積分下端の依存性として特定し、手続きの被覆率を 200 realization で実測したことです。すべての数値は記録長 4096 点・合成長 65536 点・当てはめ帯域 \(k=32\sim512\)・tol \(=0.15\) に依存します。07 節のとおり帯域を変えれば被覆率は 75.5〜86.5 % の範囲で動きます。04 節の「予言の 0.79〜0.88 倍」は \(\beta=3\)・200 realization での値です。05 節の「乗ってはいけない最小値」は本稿の合成条件でのもので、一般にどの \(\beta\) でどの \(\alpha\) が偶然当たるかは条件によって変わります。本稿も合成データのみで、実際の測定データには当てていません。 ── 印刷する場合はブラウザの「印刷」から「PDF に保存」を。

印刷 / PDF 化:⌘+P(Windows は Ctrl+P)。「答えを見る」で解答が開きます。