考える波 第 59 回 / 第 XIX 部・道具を持って 物理へ戻る
系統誤差を「見積もる」から「測る」へ
系統誤差を測る ──
そしていちばん効いたのは、窓ではなかった
第 58 回の「窓を変えると値が動く」を逆手に取り、動かない区間を探す手続きを作りました。
★★★ 途中で分かったこと ── あの「行きすぎ」は雑音ゼロでも出ます。主役は記録の分解能より低い周波数、つまり記録の中では傾きにしか見えない成分でした。
★★ だから直線を 1 本引くだけで消えます(\(\beta=5\) で \(5.512\to4.992\))── 第 58 回で「窓の階数を上げろ」と書きましたが、もっと安い手がありました。
第 58 回で符号の決まった系統誤差が出ました。★ 今回はそれを測定の道具に変えます ── 動かし方から、真の値と誤差の両方を取り出す。
01\(\hat\beta\) を \(\alpha\) について走らせる
窓の階数 \(\alpha=1\)(矩形)から 8 まで。32 realization の平均です。
| 真の \(\beta\) | \(a{=}1\) | 2 | 3 | 4 | 5 | 6 | 7 | 8 |
|---|---|---|---|---|---|---|---|---|
| 1 | 0.990 | 0.989 | 0.991 | 0.988 | 0.988 | 0.987 | 0.988 | 0.990 |
| 3 | 2.033 | 3.160 | 2.992 | 2.990 | 2.989 | 2.988 | 2.989 | 2.990 |
| 5 | 1.984 | 4.020 | 5.512 | 5.061 | 4.990 | 4.990 | 4.990 | 4.991 |
| 6 | 1.981 | 4.021 | 6.000 | 7.623 | 5.997 | 5.992 | 5.990 | 5.992 |
| 7 | 1.981 | 4.021 | 6.004 | 8.001 | 7.633 | 7.182 | 6.992 | 6.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\) | 2 | 3 | 4 | 5 | 6 | 7 | 8 |
|---|---|---|---|---|---|---|---|---|
| 1 | 0.005 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 |
| 3 | −1.009 | +0.125 | 0.001 | 0.001 | 0.001 | 0.001 | 0.001 | 0.002 |
| 5 | −3.012 | −0.988 | +0.625 | +0.068 | 0.003 | 0.003 | 0.004 | 0.004 |
| 7 | −5.012 | −2.988 | −1.000 | +1.002 | +0.918 | +0.252 | 0.007 | 0.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.0625 | 0.25 | 0.5 | 1.0 |
|---|---|---|---|---|
| 推定 \(\beta\) | 5.625 | 5.062 | 5.006 | 5.002 |
02節の中心 ── 漏れの主役が分かった
★★★ 記録の分解能 1 bin より下(0.0625〜1 bin)を入れたときだけ、行きすぎが出ます。
★★ 漏れの主役は「記録の中では傾きにしか見えない成分」でした ── 周期が記録より長い波。第 5 回・第 56 回で「観測時間より長い周期は定数と区別がつかない」と書いた、まさにその成分です。
★ 区別がつかないだけでなく、測定を汚していた。
03★★★ ならば、傾きを引けばよい
記録の中で傾きにしか見えないなら、当てはめて引いてしまえば消えるはずです。ハン窓(\(\alpha=3\))で、\(q\) 次多項式を引いてから測ります。
| 真の \(\beta\) | 引かない | \(q{=}0\)(平均) | \(q{=}1\)(直線) | \(q{=}2\) | \(q{=}3\) |
|---|---|---|---|---|---|
| 3 | 2.992 | 2.992 | 2.991 | 2.991 | 2.991 |
| 4 | 4.000 | 4.000 | 3.992 | 3.992 | 3.991 |
| 5 | 5.512 | 5.512 | 4.992 | 4.993 | 4.992 |
| 6 | 6.000 | 6.000 | 6.035 | 5.994 | 5.992 |
| 7 | 6.004 | 6.004 | 6.504 | 6.039 | 6.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 回の予言 | 比 |
|---|---|---|---|---|
| 1 | 1.000 | 0.1170 | 0.0865 | (偏りあり) |
| 2 | 1.234 | 0.1815 | 0.0961 | (偏りあり) |
| 3 | 1.500 | 0.0839 | 0.1059 | 0.792 |
| 5 | 1.944 | 0.1011 | 0.1206 | 0.838 |
| 8 | 2.473 | 0.1193 | 0.1360 | 0.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\) | 2 | 3 | 4 | 5 | 6 | 7 | 8 |
|---|---|---|---|---|---|---|---|---|
| 3 | 1.011 | 0.242 | 0.085 | 0.093 | 0.100 | 0.106 | 0.112 | 0.116 |
| 4 | 2.019 | 0.026 | 0.085 | 0.093 | 0.099 | 0.105 | 0.111 | 0.116 |
| 6 | 4.019 | 1.978 | 0.008 | 1.686 | 0.097 | 0.104 | 0.110 | 0.115 |
| 7 | 5.019 | 2.979 | 0.998 | 1.001 | 0.885 | 0.254 | 0.110 | 0.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 | 被覆率 |
|---|---|---|---|---|---|---|---|
| 1 | 1.015 | 0.129 | 0.057 | 0.078 | 0.178 | 0.317 | 87.0 % |
| 3 | 3.021 | 0.129 | 0.060 | 0.092 | 0.187 | 0.366 | 84.0 % |
| 5 | 5.014 | 0.132 | 0.062 | 0.093 | 0.201 | 0.224 | 84.5 % |
| 7 | 7.009 | 0.136 | 0.068 | 0.105 | 0.216 | 0.121 | 79.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\) | 3 | 3.021 | 0.129 | 84.0 % |
| \(k=32\ldots512\) | 6 | 6.006 | 0.133 | 83.0 % |
| \(k=64\ldots256\) | 3 | 3.027 | 0.352 | 75.5 % |
| \(k=64\ldots256\) | 6 | 6.015 | 0.367 | 78.5 % |
| \(k=16\ldots1024\) | 3 | 3.019 | 0.074 | 86.5 % |
| \(k=16\ldots1024\) | 6 | 6.008 | 0.075 | 85.0 % |
★ 帯域を変えると被覆率も動きます(75.5〜86.5 %)── 推定値そのものはほとんど動かない(3.02 前後)のに。
★★ つまりこの手続きも「仮定つき」(第 55 回)。仮定は三つ ── ① 当てはめ帯域 ② tol ③ 窓の族。手続きを作ったからといって、仮定が消えたわけではありません。
練習問題
- 02 節で「漏れの主役は記録より長い周期」と分かった。第 56 回とどうつながるか。
答えを見る
第 56 回は「見えない」と書きましたが、今回は「見えないのに効いている」が分かりました。
★ 第 56 回:「観測時間より長い \(\tau\) は折れ曲がりが見えない」。★★ 本回:それでも裾は窓を通って高周波に漏れ、傾きを \(+0.6\) も汚す。
★★★ 第 57 回 07 節の練習問題で「\(n\) が小さいほど端の外の情報が遠くまで効く」と書いたのは、これのことでした ── あのときは「原理的には \(\tau_{\max}\) を推定できるかもしれない」と書きましたが、実際に起きるのはまず汚染です。
★ 順序としては:まず引いて(03 節)、それから測る。推定に使うのはその後の話。 - 03 節で「平均を引くだけでは効かない」と出た。なぜ 1 次から効くのか。
答えを見る
窓が 0 次の漏れをすでに消しているからです。ハン窓(\(\alpha=3\))は端で値も傾きも 0 なので、一定成分(\(f=0\))の漏れは \(f^{-3}\) で十分に落ちています。
★ 残るのは \(f\) が 0 でない低周波 ── 記録の中では「直線」に見える成分。★★ 直線を引くと、その成分の大半が消えます。
★★★ 一般化すると:\(q\) 次を引くのは「記録の中で \(q\) 次多項式に見える低周波」を消すことで、それは実効的に窓の階数を上げるのと同じ働きをします。★ 違いは代金:窓の階数を上げるとばらつきが増えますが(04 節)、トレンド除去は自由度を \(q+1\) 個 減らすだけ ── 4096 点なら ほぼ無料です。 - 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 % と動きます。手続きを作っても、仮定は消えません。
★ 次回はその成分そのものを見ます ── トレンドと信号は、どこで分かれるのか。
★★ 第 5 回の「定数と区別がつかない」、第 29 回の繰り込み、そして本回のトレンド除去 ── どれも「長い方を切る」操作です。★ 切ったものはどこへ行くのかを調べます。
この文書は「考える波」シリーズ第 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 に保存」を。