考える波 第 51 回 / 第 XVIII 部・誤差棒をつける

第 49 回の手引きは、雑音のない関数でしか試していなかった

桁数に 誤差棒をつける ──
周期図 1 本の ばらつきは、いつでも 5.57 dB 第 50 回で「次に手をつけるならこれ」と書いた一行です ── 雑音のある実データで「何桁 の冪則か」を言う。
★ 周期図 1 本の ばらつきは 5.57 dB信号の強さにも 周波数にも 測定器にもよりません。
★★ そして分かったこと ── 第 49 回の手引きは、このままでは実データに使えません。1 万本 平均しても、6 桁 の冪則を 2.6 桁 としか報告できませんでした。

必要な道具:第 9, 22, 27, 42, 47, 48, 49 回計算:kenshou/gosa.js

第 27 回からずっと「何桁 続くか」で語ってきました。
★ その桁数自体に誤差をつけるのが筋です ── 第 50 回 05 節の一行めから、ここで手をつけます。

01周期図 1 本の ばらつきは、いつでも同じ

周期図の 1 本は「真のスペクトル × \(\chi^2_2/2\)」 $$\ln\!\left(\frac{\chi^2_2}{2}\right):\quad \text{平均}=-\gamma=-0.5772,\qquad \text{分散}=\frac{\pi^2}{6},\quad \text{標準偏差}=\frac{\pi}{\sqrt6}=1.2825$$
モンテカルロ理論
平均 \(E[\ln(\chi^2_2/2)]\)−0.575227−0.577216
標準偏差1.2790371.282550
★ 同じものを dB で5.5548 dB★ 5.5700 dB

01節の中心

周期図 1 本の ばらつきは 5.57 dB。
信号の強さにも 周波数にも 測定器にも よりません。

これは「測定が下手だから」ではありません ── 1 本の周期図は、実部と虚部の 2 自由度しか持たないからです。
★ 平均が \(-0.5772\)(= \(-2.51\) dB)ずれるのも有名ですが、これは切片に効くだけで傾きには効きません

02傾きの誤差は \(1.2825/\sqrt M\)

周波数が \(0\sim F\) に一様なら \(\mathrm{Var}(\ln f)=1\) ちょうど $$\ln S=c-\beta\ln f+\varepsilon\quad\Longrightarrow\quad \mathrm{sd}(\beta)=\frac{1.2825}{\sqrt M}$$
\(M\)(周期図の本数)\(\mathrm{sd}(\beta)\)(実測)\(1.2825/\sqrt M\)
2000.0895320.0906900.9872
10000.0394020.0405580.9715
200000.0092590.0090691.0209

\(\mathrm{Var}(\ln f)=1\) は「線形に並んだ周波数」の性質です ── 桁数を広げても 1 のまま。効くのは本数 \(M\) だけ。

03★ 第 42・47 回の「ずれ」を測るのに、何本 要るか

\(\text{dB/oct}=3\beta\) なので \(\mathrm{sd}(\text{dB/oct})=3.85/\sqrt M\)。

測りたい差\(1\sigma\) に要る \(M\)\(3\sigma\) に要る \(M\)試料長 \(N\)
★ 乱流の間欠性(第 42 回)0.08122.3e+32.0e+44.0e+4
★ 宇宙の \(n_s\)(第 47 回)0.10531.3e+31.2e+42.4e+4
1 階ぶんの差(3 dB/oct)1.61530

03節の中心

0.081 dB/oct を \(3\sigma\) で言うには、周期図 2 万本 = 4 万点 の試料。意外に少ない ── 統計だけなら手が届きます。
★★★ ところが実際に難しいのは、そこではありませんでした。

04★★ 系統誤差 ── まず「打ち消しあう」ことに気づく

真のスペクトルは無限に続く冪則ではなく、両端に曲がりがあります(慣性領域の幅を \(R\) 桁 とします)。その真ん中に幅 \(D\) 桁 の窓を置いて直線を当てます。

(a) 対称な模型\(D=0.5\)\(D=1\)\(D=2\)\(D=3\)
\(R=3\) 桁1.3e−143.3e−152.4e−148.0e−15
\(R=6\) 桁6.7e−161.8e−142.1e−141.1e−14

04節 (a) ── ★★★ 予想していなかった結果

ほぼ完全に 0 ── 両端の曲がりが、傾きに逆向きに効いて打ち消しあいます。

「窓を慣性領域の真ん中に置く」だけで、主要な系統誤差は消えます。
★★ これは実務でいちばん効く注意でした。本稿を書くまで気づいていませんでした。

ただし対称性に依存しています。実際の乱流は高域が指数的に切れ(散逸域)、低域は緩やかな曲がりです ── 非対称にすると打ち消しが効きません。

(b) 非対称な模型\(D=0.5\)\(D=1\)\(D=2\)\(D=3\)
\(R=3\) 桁0.09120.09820.12340.1380
\(R=4\) 桁0.03030.03320.04590.0692
\(R=6\) 桁0.00310.00340.00490.0083
(c) 窓を真ん中からずらす(対称模型、\(R=4,D=2\))★ 系統誤差[dB/oct]
ずれ 0 桁−0.00000
0.25 桁0.00937
0.50 桁★ 0.03201
1.00 桁0.26611

04節の中心

★★★ 0.5 桁 ずらすだけで 0.03 dB/oct 級 ── 打ち消しはすぐ壊れます。
慣性領域が 3 桁 しかなく非対称なら、系統誤差は 0.09〜0.14 dB/oct ── 第 42 回の 0.081 dB/oct と同じかそれ以上です。
★★ 統計誤差は \(M\) を増やせば減りますが、系統誤差は減りません。

05★★★ 窓の幅には、最適値がある

統計は窓を広げると減り、系統は広げると増える $$\text{統計}=\frac{3\times1.2825}{\sqrt{nD}\cdot\frac{\ln10}{\sqrt{12}}\cdot D}\ \propto D^{-3/2},\qquad \text{系統}=\text{04 節の表}\ \nearrow$$
慣性領域 \(R\)★ 最適な窓 \(D\)統計系統★ 合計[dB/oct]
3 桁2.95 桁0.11420.1382★ 0.1793
4 桁3.95 桁0.07370.08880.1155
6 桁5.35 桁0.04680.0355★ 0.0587

05節の中心

★★★ 最適点では、統計誤差と系統誤差がほぼ同じ大きさになります。片方だけ小さくしても総和は減りません ── だから「もっと長く測る」だけでは頭打ち。
★★ 慣性領域が 3 桁 しかないと、最良でも 0.18 dB/oct。第 42 回の 0.081 dB/oct には届きません。6 桁 あれば 0.059 まで落ちて、ようやく届きます。
つまり「長く測る」より「広い慣性領域を持つ系を選ぶ」ほうが効く ── 実験計画としては、そちらが先です。

06★ 第 49 回の手引きは、このままでは使えない

第 49 回の手順(局所的な階数が \(\pm0.1\) の最長区間)を、雑音つきでやってみます。真の信号は 6 桁 にわたる 1/f です。

平均本数 \(K\)1 点あたりの sd[dB]★ 測れた「平坦な桁数」真の値
15.5700.049 ± 0.051 桁6.0 桁
1000.5570.205 ± 0.076 桁6.0 桁
100000.056★ 2.592 ± 1.039 桁6.0 桁

06節の中心 ── ★★★ 本回の要

1 万本 平均して 1 点あたり 0.056 dB にしても、
6 桁 の冪則を 2.6 桁 としか報告できません。

原因は「隣同士の差で階数を作っていること」です。差を取ると雑音が \(\sqrt2\) 倍 になり、さらに間隔で割るので \(1/\Delta\) 倍 に増幅される。
★★ 第 49 回の手引きは、雑音のない関数でしか試していませんでした ── 実データに当てれば破綻します。
これは本シリーズが自分で見つけた欠陥です。第 50 回で「次に手をつける」と書いた通りにやったら、出てきました。

07★ 直し方 ── 差ではなく、窓ごとの当てはめで測る

窓の幅 \(w\)[桁]階数の sd★ 測れた「平坦な桁数」上限 \(6-w\)
0.100.698410.051 ± 0.017 桁5.9
0.300.127950.302 ± 0.087 桁5.7
0.600.044663.409 ± 1.328 桁5.4
1.000.02005★ 5.017 桁5.0
1.500.011284.535 桁4.5

07節の中心

局所的な階数を「隣との差」でなく「幅 \(w\) の窓での最小二乗」で作ると、直ります。窓の中の点数だけ雑音が \(\sqrt{\ }\) で減るから。
★★★ ただし幅 \(w\) の窓は両端で \(w/2\) ずつ使えなくなるので、報告できる桁数の上限が \(6-w\) に下がります。実際 \(w=1.0\) で 5.017(上限 5.0)、\(w=1.5\) で 4.535(上限 4.5)── 端で頭打ちです。
★★ 狭すぎると雑音で途切れ、広すぎると端を失う ── 窓の幅そのものが、第 9・22 回の「分解能の配分」でした。

08第 49 回の手引きへの追記

手順注意
③-a 1 点の誤差を 5.57 dB と置く周期図ならいつでもこの値
③-b \(K\) 本 平均したら \(5.57/\sqrt K\) dB平均は周波数方向でも試行方向でもよい
★ ③-c 局所的な階数は「窓での当てはめ」で★ 隣との差は使わない(06 節)
★ ③-d 窓の幅 \(w\) は 0.5〜1 桁★ 報告できる桁数は \(D_{\rm total}-w\) が上限
★ ③-e 窓は慣性領域の真ん中に置く★ 端の曲がりが打ち消しあう(04 節)
③-f 統計と系統が釣り合う \(M\) で止めるそれ以上 測っても精度は上がらない

09波として言えること・言えないこと

項目判定根拠
★ 周期図 1 本の ばらつき◎ 数値★ 5.57 dB(何によらず一定)
\(\mathrm{sd}(\beta)=1.2825/\sqrt M\)◎ 数値★ 桁数ではなく本数で決まる
0.081 dB/oct に要る試料長◎ 数値\(3\sigma\) で 4 万点(統計のみ)
★ 対称な窓では系統誤差が打ち消す◎ 数値★ 真ん中に置くだけで消える
0.5 桁 ずらすと壊れる◎ 数値0.032 dB/oct
★★ 窓の幅に最適値がある◎ 数値★ そこで統計と系統が同じ大きさ
★ 素朴な桁数推定は破綻する◎ 数値★ 1 万本 平均しても 2.6 桁
窓での当てはめで直る◎ 数値\(w=1\) で 5.0 桁(上限どおり)
★ 実データでの検証× 未着手★ 本稿も合成データのみ

練習問題

  1. 01 節の「5.57 dB」は、なぜ測定器によらないのか。
    答えを見る
    周期図の 1 本が、実部と虚部の 2 自由度しか持たないからです。どんなに良い測定器でも、ある周波数の 1 回の推定値は複素数 1 個ぶんの情報しかありません。
    ★ その大きさの二乗は \(\chi^2_2\) に従い、その対数の分散は \(\pi^2/6\) という普遍定数
    ★★ 減らす方法は一つだけ ── 独立な推定値を増やして平均すること(\(K\) 本 で \(5.57/\sqrt K\) dB)。周波数方向に束ねても、試行を繰り返しても同じです。★ 第 9 回の「配分は選べるが総量は動かせない」と同じ構図で、ここでは「自由度の数」が総量にあたります。
  2. 04 節で系統誤差が打ち消しあった。実データでこれを使うには何に気をつけるか。
    答えを見る
    三つあります。
    ★ ① 両端の曲がりが同じ形か:04 節 (b) のとおり、片方が指数的に切れると打ち消しは効きません。乱流の散逸域はまさにそれです。
    ★ ② 窓が本当に真ん中か:0.5 桁 ずれると 0.03 dB/oct(04 節 (c))。慣性領域の端の位置を独立に知っている必要があります。
    ★ ③ 打ち消しを「当てにしない」こと打ち消しは幸運であって保証ではないので、ずらしたときにどれだけ動くかを必ず調べる(感度解析)。
    ★★ 本稿の立場では、打ち消しは「うまくいけば得をする」であって「前提にしてよい」ではありません。
  3. 06 節で手引きの欠陥が見つかった。これは第 49 回が間違いだったということか。
    答えを見る
    間違いではなく、適用範囲の宣言が足りていました。第 49 回の本文にも「本稿は雑音のない関数だけを扱っており、雑音のある実データでは桁数の推定に統計誤差が入る」と書いてあります。
    ところが実際にやってみると、入るのは「統計誤差」ではなく「大きな系統的な過小評価」でした ── 誤差棒がつくのではなく、答えそのものが 6 桁 から 2.6 桁 へ落ちる。
    ★★ これは第 49 回の但し書きが甘かった、ということです。「誤差が入る」と「使えない」は別で、実際は後者でした。
    第 15 回から「負けた計算を消さない」と決めている以上、こう書くのが筋です ── そして直し方(07 節)も一緒に出せたので、差し引きでは前進です。

まとめ 誤差棒は、二つに分かれていた

周期図 1 本の ばらつきは 5.57 dB(\(\pi/\sqrt6=1.2825\) nats)── 信号の強さにも周波数にも測定器にもよりません。平均の \(-2.51\) dB のずれは切片に効くだけで、傾きには効きません。

傾きの誤差は \(\mathrm{sd}(\beta)=1.2825/\sqrt M\)(数値で確認)── 桁数ではなく本数で決まります。これだけなら、0.081 dB/oct を \(3\sigma\) で言うのに 4 万点。統計だけなら手が届きます。

★★★ ところが ── 窓を慣性領域の真ん中に置くと、
両端の曲がりが逆向きに効いて系統誤差が打ち消しあいます(\(10^{-14}\) まで)。
★ そして 0.5 桁 ずらすだけで 0.032 dB/oct に戻ってしまう。

非対称な場合(高域が指数的に切れる = 実際の乱流)は打ち消しが効かず、慣性領域 3 桁 では 0.09〜0.14 dB/oct。統計と系統を合わせると窓の幅に最適値があり、最適点では両者がほぼ同じ大きさになります ── \(R=3\) 桁 で最良 0.18、\(R=6\) 桁 で 0.059 dB/oct。「長く測る」より「広い慣性領域を持つ系を選ぶ」ほうが効きます。

★★★ そして本回の要 ──
第 49 回の手引きは、このままでは実データに使えませんでした。
1 万本 平均して 1 点 0.056 dB にしても、6 桁 の冪則を 2.6 桁 としか報告できない。

原因は「隣同士の差で階数を作っていたこと」── 差は雑音を \(\sqrt2/\Delta\) 倍 に増幅します。直し方は「幅 \(w\) の窓での最小二乗」で、\(w=1\) 桁 なら 5.0 桁 まで復元できました。
★★ ただし幅 \(w\) の窓は両端を \(w/2\) ずつ食うので、報告できる桁数の上限が \(D_{\rm total}-w\) に下がります ── 狭すぎれば雑音で途切れ、広すぎれば端を失う。窓の幅そのものが、第 9・22 回の「分解能の配分」でした。

正直に本稿も合成データだけで、実際の測定データには当てていません。★ そして第 49 回の但し書きは甘かった ── 「誤差が入る」ではなく「使えない」でした。

この文書は「考える波」シリーズ第 51 回です。計算は kenshou/gosa.js で行っています。本回の数値について:対数周期図の統計量のモンテカルロ、傾きの標準偏差、必要な試料長、系統誤差の計算、桁数推定のシミュレーション、窓つき当てはめによる改善は、すべて本稿で計算したものです。周期図が \(\chi^2_2\) 分布に従うこと、対数周期図の平均が \(-\gamma\)・分散が \(\pi^2/6\) であること(ゲウェケ=ポーター=ハダック法の基礎)、およびそれによる傾き推定の誤差 ── これらはいずれもよく知られた結果で本稿の発見ではありません ── 本稿がしたのは、それを本シリーズの「dB/oct」と「桁数」という言葉に翻訳し、第 49 回の手引きに当てて破綻を見つけ、直したことです。04 節の「系統誤差が打ち消しあう」は、両端の曲がりが同じ関数形で窓が正確に中心にある場合に限った結果であり、一般には成り立ちません((b)(c) がその反例です)。05 節の最適化は「1 桁あたり \(n=100\) 本の独立な推定値」という理想化のもとでの計算で、実際には対数ビン化により低周波側の本数が少なくなります。03 節の試料長は白色に近い過程を仮定した目安で、実データでは窓関数・漏れ・前処理の影響が加わります。本稿も合成データのみで、実際の測定データには当てていません。 ── 印刷する場合はブラウザの「印刷」から「PDF に保存」を。

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