前置き

私はほぼ毎日自転車（ママチャリ）に乗って移動する。日頃通う道はそれほど整備されておらず、自転車レーンはないため道路脇を走る。交通量はそれほど多くはないが、一分も走れば何台かの車には追い越される。¹ 私は割と心配症なので、事故に遭うのを避けるため、道路脇を走行するときは定期的に振り返り、後方確認をする。これは私が個人的に「車による自転車追い越しのプロトコル」と呼ぶ、以下の手順を踏んで追い越してもらうのが最も安全だと考えているからだ。

車による自転車追い越しのプロトコル

自転車側: 後方を確認し車を視認する。

車側: 自転車が後方確認したことを視認する。

車が追いつくと推測されるタイミングの少し前に、必要があれば自転車が少し脇に寄るか、車が対向車線に少しはみ出す。

車が自転車を追い越す。

このように運転手との間で緩やかな相互認証を行うことにより、単に追い越してもらうよりも安全になると考えている。後方確認の別のメリットとして、車が来ないことが保証された安全な時間が確保できる点がある。つまり、一度振り向いて車がいなければ、しばらくは振り向かずとも車がすぐそばにいないことがわかる。この間は道路の中心寄りを走っても安全なので、快適に走行でき、やや危険性の高い路肩（例えば、ガタガタしていたり、伸びた雑草が道路に飛び出していたりする）から距離を取ることができる。

ただし、後方確認には以下のようなデメリットもある。

短時間に後方を向き、前方を向き直す必要があるため、首をある程度勢いよく捻る必要があり、怪我の危険性がある。特に寝違えているときは一層辛い。
振り向きの動作は上体の旋回を伴い、自転車の進行方向が旋回方向にずれ、車道の中央側に飛び出してしまう。
前方への注意が途切れる。

この二、三点目は特に重要で、後方確認自体が交通事故の可能性を高めてしまう可能性がある。つまり、後方確認は常に有益とは限らない。大雑把に言って、車が後方から迫ってきたとき、後方確認をする、しないによって、以下のような4つの結果があり得るだろう。

後方確認をしない
1. 車がうまく回避して追い越す
2. 車側が自転車が急に飛び出さないか不安になり、不規則な事態へ
後方確認をする
1. プロトコルが機能し、車が安全に追い越す
2. 車が直後にいるタイミングで後方確認し、自転車が道路中央側に飛び出し不規則な事態へ

どちらの場合にも良い結果と悪い結果があり、どちらが起こるかには不確実性がある。とはいえ、状況に応じてより適した選択はあるだろう。すなわち、後方確認をしない方が良い結果が起こる可能性が高い状況、またはその逆もある。

例えば、

自転車レーンが整備され、十分なスペースがある場合、後方確認をせず、まっすぐに進むのが安全
交通量が多い場合、後方確認に伴う自転車の道路中央側への飛び出しの危険性が高まるので、確認しないのが安全
交通量がそれほどでもない場合、後方確認自体の危険性は低く、確認することで上記のプロトコルが機能し、安全性が増す

後方確認が不適切な状況においては私たちができることはなく、単に真っ直ぐに走行し、事故に遭わないことを祈るのみだ。² よって、この文章では、それほど交通量が多くない道において、どの程度後方確認をするのがいいかを考える。

適切な後方確認頻度とは

まず、後方確認は可能な限り少ない回数にするべきである。なぜなら上記のように、首への負荷や飛び出しの危険性があるからだ。一方で確認を怠りすぎると、確認をしない場合のデメリットを受ける可能性が高まるし、不適切なタイミングでの確認による飛び出しの危険性がある。よって、安全を保証しつつも可能な限り回数を抑えるのがいい。それはどれくらいの頻度だろうか？

これは細かいことを無視すれば、 $T-t$ で計算するのがいいだろう。ここで、

$T$ : 次に車が自分に追いつきうる最短予想時間
$t$ : 後方の状況に対応するために必要な猶予時間

である。

より詳しくいうと、 $T$ は、直前の確認で得た後方の情報、その道路を走る車の典型的な速度、自転車の速度などから推測される、車が次に自分に追いつきうるまでの最短時間である。

例えば、最も単純な場合として、脇道のない一直線の道を走る場合、直前に振り返ったときに見えた車までの距離を $d$ としよう。車がいなかった場合には視認できる最大距離を $d$ とする。そして、その道を走る車の（予想）速度を $v_{C}$ 、自転車の速度を $v_{B}$ としたとき、 $T = d / (v_{C} - v_{B})$ で求められる。

猶予時間 $t$ は後方確認し、対応を取るのに必要な時間である。例えば、

t = (後方確認に要する時間) + (後方から車が接近していた場合、脇に寄るのに必要な時間)。

この猶予時間 $t$ も一般には車と自転車の速度等に応じて変化するだろうが、それほど大きくは変わらないだろう。よって今回の考察では $t$ は定数と考える。私の個人的な体感では $t$ は4秒くらいと考えれば十分だと思う。

この時間 $T-t$ 以内に確認を繰り返すことで、常に後方の状況を対処可能な状態に保つことができる。もちろんイレギュラーな事態が起こらない限りではあるが。

猶予時間 $t$ は定数としたので、 $T$ がわかれば、適切な後方確認頻度が求められる。上では単純な場合の $T$ の計算例を与えたが、一般には通る道の特性に応じた考察が必要だろう。

例えば、複数の脇道からの流入がある道路を走る場合、先ほど求めたように直進して追いつく車だけでなく、脇道から入ってきて追いついてくる車の存在も考慮する必要がある。この場合、先ほどのように $T = d/(v_{C} - v_{B})$ とするのではなく、脇道から車が入ってきた場合に、それが自分に追いつくまでの（予想）時間を $T'$ として、 $d/(v_{C} - v_{B})$ と $T'$ の小さい方を $T$ として採用する必要があるだろう。また脇道から入ってきた車は減速しているため、加速して $v_{C}$ に達することを考慮して $T'$ を求める必要がある。このため、暗算できるような簡単な公式という形では結論を出すことはできないが、大して精密さを要求される計算でもないから、例えば等加速度運動と考えれば、高校物理で扱うような単純な計算で求めることができるだろう。

他にも状況に応じて、 $T$ の定め方は多様にあるだろうが、ここでは典型的な場合の一つとして、今述べたような、直線と脇道から車が追いついてくる場合のみを考える。

具体例における計算

具体例で計算を行う。付録として、この計算を行うコードを文章の末尾に載せた。実際に私のいつも通っている道をモデル化して考える。

私が通る道はそれほど交通量が多くなく、それほど広くもない。おそらく車の時速は40km/hぐらいで、私の走行速度は15km/hぐらい。見通しはそれなりによく、200mぐらいは見通せるが、いくつか脇道が合流しているので、そこから車が入ってくる可能性を考えると、一回の振り向きで保証できる安全な距離の最大は大体80mぐらいだろう。脇道から合流する車は10km/hまで速度が落ちているとし、40km/hに達するまでに達するまでの加速度0.2g (gは重力加速度)とする。³ まず、直進してくる車が追いつく時間を考える。直前に振り向いた際には後方に車がいなかったとすると $d$ は200mなので、 $d/(v_{C} - v_{B}) = 200/7$ ≒ 29秒。

次に脇道から入ってきた車が追いつく時間を考える。簡単な計算によって、脇道から入ってきた10km/hの車が40km/hに達するまでの時間は大体4秒でそれまでの走行距離は18mとわかる。時速を秒速に直すと、15km/h ≒ 4m/s. 40km/h ≒ 11m/h。相対速度は25km/hで秒速に直すと約7m/s。よって、80m 後方の脇道から入ってきた車が追いつくまでの時間は、 $T' = (80 + 4 \cdot 4 - 18)/7 + 4$ ≒ 15秒。

$T$ は29と15の小さい方なので15秒。猶予時間 $t$ は4秒とすると、結局 $T - t$ は11秒。よって、11秒に一回振り向くというのが良さそうである。ただしこれは猶予時間ギリギリの振り向きになるので、もう少し余裕を持って9秒に1回ぐらいがちょうどいいかもしれない。

実践と反省

この計算結果を踏まえ、今日実際に、脳内で十秒数えるごとに振り返りながら走行してみた。感想としては、この振り返り頻度は想像よりもずっと頻繁で煩わしく、やや異常な振る舞いだと感じた。もし、このようにしきりに振り返る他人を見たら、私は何らかの神経症的な性向を疑う。

というわけで、早々に振り返り頻度を下げ、そのまましばらく走っていると、後ろからエンジン音が聞こえ、振り返るとそれなりに近く（十メートル程度）まで車が迫っていた。これは余裕をもって対処はできたのだが、より早く確認していれば一層安全だったわけで、十秒に一回は多すぎるにしても、秒数を決めての規則的な振り返り自体は真剣に考えるべきかもしれない。

付録: コード

計算に使ったコードを載せておく。「ユーザー入力」の部分の値を変更することで、さまざまな状況下に対応できる。そう難しい内容でもないが、コーディングはChatGPTに頼った。一通り目を通して正しい動作をすることは確認したが確実ではない。

このコードを通じて、sympy というのを初めてみた。これは不定元を代数的に扱う機能を備えたライブラリのようで、例えば、二次方程式の係数を未知定数にしたままで、解の公式を用いて解の表示を与えることができ、後から係数に値を代入することができる。便利なものだと感じた。

コードを開く

from sympy import symbols, Eq, solve

def kmph_to_mps(speed):
    return speed * 1000 / 3600

# 定数と変数の定義
t = symbols('t')
V_b, v_l, v_e, a, d = symbols('V_b v_l v_e a d')
g = 9.8

# ユーザー入力
V_b_input = 15  # 自転車の速度(km/h)
v_l_input = 10  # 流入する車の初速度(km/h)
v_e_input = 40  # 車の終端速度(km/h)
a_input = 0.2   # 車の加速度(G). 重力加速度を単位として用いている
d1_input = 80   # 脇道から車の現れうる最短距離(m)
d2_input = 200  # 直進してくる車との最短距離, または後方の最大見通し(m)

# 単位変換
V_b_val, v_l_val, v_e_val = map(kmph_to_mps, [V_b_input, v_l_input, v_e_input])
a_val = a_input * g
print(f"V_b:{V_b_val:.3f} m/s, v_l:{v_l_val:.3f} m/s, v_e:{v_e_val:.3f} m/s, a:{a_val:.3f} m/s^2")

# 自転車と自動車の位置関数
x_bike = V_b * t
x_car_accel = -d + v_l * t + 0.5 * a * t**2
t_end_accel = (v_e - v_l) / a
x_car_const = x_car_accel.subs(t, t_end_accel) + v_e * (t - t_end_accel)

# それぞれの追いつく時刻
catch_up_time_accel = solve(Eq(x_bike, x_car_accel), t)[1]
catch_up_time_const = solve(Eq(x_bike, x_car_const), t)[0]

# 代入して計算
values = {V_b: V_b_val, v_l: v_l_val, v_e: v_e_val, a: a_val, d: d1_input}
actual_time_accel = catch_up_time_accel.subs(values)
actual_time_const = catch_up_time_const.subs(values)
actual_time_straight = d2_input / (v_e_val - V_b_val)

# 出力
if actual_time_accel < t_end_accel.subs(values):
    if actual_time_accel < actual_time_straight:
        print("脇道からの車が加速中に追いつく.")
        T = actual_time_accel
    else:
        print("直進する車が追いつく.")
        T = actual_time_straight
elif actual_time_accel < actual_time_straight:
    print("脇道からの車が等速に達してから追いつく.")
    T = actual_time_const
else:
    print("直進する車が追いつく.")
    T = actual_time_straight
    
print(f"追いつく時刻は{T:.3f}s.")

V_b:4.167 m/s, v_l:2.778 m/s, v_e:11.111 m/s, a:1.960 m/s^2 脇道からの車が等速に達してから追いつく. 追いつく時刻は14.071s.

「車」ではなく「自動車」の方が正確だが、「自転車」と字面が見分けにくいので、この文章では「車」を用いる。↩
私はそのような道を自転車で通ることは、日常におけるリスクとしては大きすぎるのではないかと思っている。生活の都合上避けられない場合もあるけれども。↩
加速度については次のページを参考にした: 車が静止状態から発進するときの加速度を教えてくださいませんか？現実の状況、レースではありません。 - Quora ↩

2023-09-11

分散分析(ANOVA)のF検定の導出について２——尤度比検定との等価性

前置き
準備
ガウス・マルコフモデル
尤度比検定とF検定の等価性

前置き

前回の文章では, 分散分析の等平均性の検定の検定統計量がF分布に従うことを示した.

kihanto.com

そこでは検定統計量は最小二乗値を用いて表された. 今回の文章の目的は, 尤度比検定を出発点にして検定を構成しても, 同じ最小二乗値を用いた検定統計量が自然に現れることを示すことである.

なぜこのようなことを考えるのか, 目的を明確に述べるために, まず前回の導出の流れを振り返る:

ガウス確率変数の二次形式はフィッシャー・コクランの定理から, 分布, 独立性が容易にわかる.
ガウス・マルコフモデルにおける, 制約条件つきと条件なしの最小二乗値( $R_{1}^{2}$ と $R_{0}^{2}$ は, ガウス確率変数の二次形式で表せる.
フィッシャー・コクランの定理から, $R_{1}^{2} - R_{0}^{2}$ と $R_{0}^{2}$ は独立な二次形式となり, 前者は非心カイ二乗分布, 後者はカイ二乗分布に従う.
制約条件が実際に成り立つとき, $R_{1}^{2} - R_{0}^{2}$ は非心度0のカイ二乗分布に従うので, 比をとることでF検定が導出できる.
分散分析モデルはガウス・マルコフモデルの特別な場合であるため, 以上の内容が適用できる.
等平均性の仮説は線形制約条件で表すことができ, その仮説に対し, $R_{1}^{2} - R_{0}^{2}$ と $R_{0}^{2}$ を具体的に計算すると, よく知られた平方和の式に行き着く.
よって, 分散分析において用いられている検定量は等平均性の仮定の下ではF分布に従う.

この導出は線形代数やガウス確率変数に関する一般的で比較的単純な結果の応用として示すもので, 論理はクリアだが, 不満点もある. それは, 最小二乗値から構成した検定統計量を考える明確な動機づけに乏しく, やや天下り的なアプローチに見えるという点である.

そこでここでは, 尤度比検定を出発点として検定を構成しても, 前回考えた最小二乗値に関する検定統計量が現れることを示す. これにより, 最小二乗値による検定統計量を考えるに至る自然な道筋が存在することを示すのが, この文章の目的である.

ただし, 統計量が「自然に導かれる」ことは, それが良い統計量だということを必ずしも意味しないことを注意しておく. 統計量の良さは検出力や分散の小ささといった目的に沿った具体的な指標で評価されるべきだからだ.

ここからの内容のアウトラインを説明する. まず, 線形代数から, 射影と一般化逆行列に関する幾つかの命題を準備する. 次に, ガウス・マルコフモデルを説明する. ここは前回の引き写しである. そして, ガウス・マルコフモデルに対し, 線形制約条件を帰無仮説とした場合の尤度比検定を考えることで, 最小二乗値を用いた検定統計量を考えることに自然に導かれることを示す.

準備

以下では, 線形代数の標準的と思われる表記法は断りなく用いる. ベクトル空間は常に適切な次元のユークリッド空間とし, 内積は常に標準内積を考えるものとする. また, 行列とそれの定める線形変換を同一視する. つまり, 行列 $A$ から定まる線形変換 $x \mapsto Ax$ のことを線形変換 $A$ と書いたり, その変換の像のことを $\mathrm{Im} \ A$ と書いたりする. 行列 $A$ の転置を $A'$ とかく.

まず前回も用いた, 線形代数から2つの命題を思い出す.

命題 1 行列 $A$ に対し, $\mathrm{Ker} \ A' = (\mathrm{Im} \ A)^{\perp}$ .

証明

ベクトル $x,y$ に対し, $x \in \mathrm{Ker} \ A'$ なら, $x'Ay = (A'x)'y = 0$ . $y$ は任意に取れるので, $x \in (\mathrm{Im} \ A)^{\perp}$ . 逆に, $x \in (\mathrm{Im} \ A)^{\perp}$ なら, $(A'x)'y = x'Ay = 0$ . $y$ は任意より, $A'x = 0$ . つまり, $x \in \mathrm{Ker} \ A'$ . (証明終わり)

線形変換 $P$ が $P^2 = P$ を満たすとき, 射影という. さらに, $P' = P$ を満たすとき, 直交射影という. 部分空間 $U$ が射影 $P$ の像であるとき, $P$ を $U$ への射影という.

命題 2 $P$ を部分空間 $U$ への直交射影とするとき, 任意の $y \in \mathbb{R}^n$ に対し,

$\min_{z \in U}|y - Pz|^2 = |y - Py|^2 = y'(I - P)y. \label{}$

証明

任意の $z \in U$ に対し,

$|y - Pz|^2 = |y - Py|^2 + |P(y-z)|^2 + 2(y-Py)'P(y-z). \label{}$

ここで, $(I-P)P = 0$ より, $(y-Py)'P(y-z) = y'(I-P)P(y-z) = 0$ . よって, $|y - Pz|^2 \geq |y - Py|^2$ . (証明終わり)

次に一般化逆行列について基本的な命題をいくつか示す. $m \times n$ 行列 $A$ に対し, $n \times m$ 行列 $B$ が $ABA = A$ を満たすとき, $B$ は $A$ の一般化逆行列という.

注意: Wikipedia に一般化逆行列の記事はあるが, やや制限の強い定義を採用しているようである. 英語版では, ここと同じ定義で扱われている.

一般化行列は常に存在する.

命題 3

任意の $m \times n$ 行列 $A$ に対し, その一般化逆行列は存在する.

証明

$\mathrm{rank} \ A = r$ とすると, 適当な正則行列 $P$ , $Q$ によって,

とできる. ただし, $I_r$ は $r \times r$ の単位行列で, それぞれの $O$ は適当なサイズの全ての成分が $0$ の行列である. $A^- := QR'P$ と定めると,

$AA^-A = (P^{-1}RQ^{-1}) (QR'P) (P^{-1}RQ^{-1}) = P^{-1}RR'RQ^{-1} = P^{-1}RQ^{-1} = A. \label{}$

よって, $A^{-}$ は一般化逆行列である. (証明終わり)

一般化逆行列は一意的とは限らない. 例えば,

とするとき, 任意の実数 $a,b,c$ に対し,

は一般化逆行列である. ただし明らかなように, $A$ が正方行列で正則なときには, 通常の逆行列が唯一の一般化逆行列である.

以下では行列 $A$ に対し, 一般化逆行列をどれでもいいので一つ固定し, $A^{-}$ と書く. 特別な一般化逆行列を取る必要があるときはその都度言及する.

命題 4 行列 $A$ に対し, 以下が成立:

$\mathrm{rank} \ A^-A = \mathrm{rank} \ A$ .
$\mathrm{Im}\ A' = \mathrm{Im} \ A'A$ .

特に, $\mathrm{rank} \ A'A = \mathrm{rank} \ A' = \mathrm{rank} \ A$ .

証明

(1): $\mathrm{rank} \ A = \mathrm{rank} \ AA^-A \leq \mathrm{rank} \ A^-A \leq \mathrm{rank} \ A$ よりわかる.

(2): ベクトル $x$ に対し, $A'Ax = 0$ なら, $0 = x A'Ax = |Ax|^2$ より, $x \in \mathrm{Ker} \ A$ . よって, $\mathrm{Ker} \ A = \mathrm{Ker} \ A'A$ なので, $\mathrm{Im} \ A' = (\mathrm{Ker} \ A)^{\perp} = (\mathrm{Ker} \ A'A)^{\perp} = \mathrm{Im} \ A'A. \label{}$
よって, 証明できた. (証明終わり)

一般化逆行列を使うと, 連立一次方程式の解をシンプルに表せる.

命題 5 $A$ を $m \times n$ 行列, $y \in \mathbb{R}^m$ とし, $y \in \mathrm{Im} \ A$ を仮定する. このとき, 連立一次方程式

$Ax = y \label{}$

の解の集合は

$\{ A^-y + (I-A^-A)u \mid u \in \mathbb{R}^n \} \label{}$

である.

証明

$y \in \mathrm{Im} \ A$ より, ある $z \in \mathbb{R}^n$ が存在し, $y = Az$ . よって, 任意の $u \in \mathbb{R}^n$ に対し,

$A(A^-y + (I-A^-A)u) = AA^-Az + (A - AA^-A)u = Az = y. \label{}$

よって, $A^-y + (I-A^-A)u$ は $Ax = y$ の解である. 特に, $u = 0$ とすれば, $A^-y$ は特殊解であるから, あとは $I-A^-A$ が $\mathrm{Ker} \ A$ への射影であることを示せばよい.

$(I-A^-A)^2 = I-A^-A$ , $\mathrm{Im} \ (I-A^-A) \subset \mathrm{Ker} \ A$ は明らかである. よって, $\mathrm{rank} \ (I-A^-A) = \mathrm{dim} \ (\mathrm{Ker} \ A)$ を示せば十分である. $I-A^-A$ は射影なので, $n = \mathrm{rank} \ A^-A + \mathrm{rank} \ (I - A^-A)$ . よって, 命題4より,

$\mathrm{rank} \ (I-A^-A) = n - \mathrm{rank} \ A^-A = n - \mathrm{rank}A = \mathrm{dim} \ (\mathrm{Ker} \ A). \label{}$

よって, 証明できた. (証明終わり)

一般化逆行列は, 射影に関するシンプルな表示も与える.

命題 6 行列 $A$ に対し, $A(A'A)^-A'$ は $\mathrm{Im} \ A$ への射影. また $(A'A)^-$ が対称なら, $A(A'A)^-A$ は直交射影.

証明

$P := A(A'A)^-A'$ とおく. このとき, $\mathrm{Im} \ A' = \mathrm{Im} \ A'A$ より, 任意の $x$ に対し, ある $z$ が存在し, $A'x = A'Az$ . これより, $P^{2}x = A(A'A)^{-}A'A(A'A)^{-}A'Az = A(A'A)^{-}A'Az = Px$ . よって, $P^2 = P$ . 後半の主張は自明. (証明終わり)

ガウス・マルコフモデル

ガウス・マルコフモデルとは, 観測値がいくつかの説明変数の線型結合と誤差との和によって与えられるとするモデルである. より具体的には, $n$ 個の観測値 $Y_1, \dots, Y_n$ が与えられたとき, それぞれの観測値は $p$ 個の説明変数 $x_{i1},\dots, x_{ip} \ (1 \leq i \leq n)$ の線形関係によって生成されると考え, そして, 観測には誤差が伴うとしたものである. また, 誤差は独立同分布で分散 $\sigma^2$ の正規分布に従うとする. つまり, 以下の観測値ベクトル $Y$ と計画行列 $X$ は, あるパラメータベクトル $\beta$ と誤差ベクトル $\varepsilon$ に対し, 以下のような関係を満たすとする:

$Y = X \beta + \varepsilon. \label{}$

ここで,

で, $(\varepsilon_i)_{1 \leq i \leq n}$ は独立同分布で平均 $0$ , 分散 $\sigma^2$ の正規分布に従うとする. $Y$ と $X$ は既知の値として扱い, パラメータ $\beta$ と分散 $\sigma^2$ は未知の値である. $r := \mathrm{rank} \ X \leq n-1$ を仮定する. これは全く制限的な仮定ではない. また, $X$ がフルランクであること, つまり, $\mathrm{rank}\ X = p$ は仮定しない.

$p \times k$ 行列 $H$ 行列は $\mathrm{rank}\ H = k$ かつ $\mathrm{Im}\ H \subset \mathrm{Im}\ X'$ を満たすとし, また $d \in \mathrm{Im}\ H' \subset \mathbb{R}^k$ とする. このとき, 以下の帰無仮説 $H_0$ , 対立仮説 $H_1$ の検定を考える:

$H_0: \ H'\beta = d, \quad H_1: \ H'\beta \neq d \label{}$

尤度比検定とF検定の等価性

帰無仮説 $H_0$ , 対立仮説 $H_1$ の尤度比検定から検定量を導出することを考える. つまり, 以下の検定統計量を考える:

$\lambda(Y) := \frac{\sup_{\beta: H'\beta = d, \sigma^2} L(\beta,\sigma^2|Y)}{\sup_{\beta, \sigma^2} L(\beta,\sigma^2|Y)}. \label{eq05}$

ここで, $L$ は尤度関数. つまり,

$L(\beta,\sigma^2|Y) := \frac{1}{(2\pi\sigma^2)^{n/2}} \mathrm{e}^{-\frac{|Y-X\beta|^2}{2\sigma^2}}. \label{}$

まず, 分母は制約条件の無い尤度の最大化であるから, 例えば, 偏微分して極値を求めることで容易に以下がわかる:

$\sup_{\beta, \sigma^2} L(\beta,\sigma^2|Y) = \left( \frac{n}{2\pi \mathrm{e}} \right)^{n/2} |Y'(I - P_X)Y|^{-n/2}. \label{}$

ここで, $P_X$ は $X$ の像への直交射影. 後の便宜のため, $Y'(I - P_X)Y = (Y - X\beta)'(I-P_X)(Y-X\beta)$ であることを注意しておく. なぜなら, $(I-P_X)X = 0$ だからである.

よって, あとは分子の値を具体的に求めれば, 尤度比が求まる. これは制約条件付きの最大化問題であるからラグランジュの未定乗数法ですぐにできそうである. しかし, ラグランジュ関数の停留点を求める方程式の解き方が分からなかった. そこで, ここではより素朴に, 制約条件を満たす $\beta$ の範囲を先に求め, その中で(制約のない)通常の最大化問題を考える. $H'\beta = d$ を満たす $\beta$ を $1$ つ固定し, $\gamma$ とする.

定理 1 以下が成り立つ:

$\sup_{\beta: H'\beta = d, \sigma^2} L(\beta,\sigma^2|Y) = \left( \frac{n}{2\pi \mathrm{e}} \right)^{n/2} |(Y-X\gamma)'(I - P_{X|H})(Y-X\gamma)|^{-n/2}$

ただし, $P_{X|H}$ は直交射影で, 像は

$\{X\alpha \mid \alpha \in \mathrm{Ker} \ H'\}. \label{eq07}$

また, $\mathrm{rank} \ P_{X|H} = r - k$ が成り立つ.

証明

$\mathrm{rank}\ H' = k$ より, $\mathrm{Ker} \ H'$ の基底 $g_1,\dots,g_{p-k}$ を $1$ つ固定し, それらを並べた行列を $G := (g_1, \dots, g_{p-k})$ とすると, 連立一次方程式 $H'\beta = d$ の任意の解は $\gamma + G\delta \ (\delta \in \mathbb{R}^{p-k})$ と書ける. よって,

$\sup_{\beta: H'\beta = d} L(\beta,\sigma^2|Y) = \sup_{\delta, \sigma^2} \frac{1}{(2\pi\sigma^2)^{n/2}} \exp \left(-\frac{|Y-X\gamma - XG \delta|^2}{2\sigma^2}\right). \label{}$

右辺は制約条件の無い関数の最大値である. 右辺の $\sup$ の中身の関数の対数を $\ell$ と書き, $\frac{\partial}{\partial \delta}\ell = 0$ とすると,

$G'X'XG \delta = G'X'(Y-X\gamma). \label{}$

命題4より, 右辺は $\mathrm{Im} \ G'X'XG$ に含まれるので, これを満たす $\delta$ は,

$(G'X'XG)^-G'X'(Y-X\gamma) + (G'X'XG)^-(G'X'XG)u \quad (u \in \mathbb{R}^{p-k})$

と書ける. ただし, $(G'X'XG)^-$ は $G'X'XG$ の一般化逆行列の $1$ つで対称なものとする. 対称なものが取れることは例えば, $G'X'XG$ のスペクトル分解を $a_1p_1 + \cdots a_mp_m \ (a_i > 0, \ p_i$ は互いに直交する直交射影) とするとき, $(G'X'XG)^- = a_1^{-1}p_1 + \cdots a_m^{-1}p_m$ ととれることからわかる. $P := XG(G'X'XG)^-G'X'$ と定めると, 命題6より, $P$ は $\mathrm{Im} \ XG$ への直交射影であり, $u \in \mathbb{R}^{p-k}$ に依らず,

$XG\delta = P(Y-X\gamma) \label{}$

である.

次に, $\frac{\partial}{\partial \sigma^2}\ell = 0$ を解くと,

$\sigma^2 = \frac{(Y-X\gamma)'(I-P)(Y-X\gamma)}n. \label{}$

よって,

$\sup_{\beta: H'\beta = d, \sigma^2} L(\beta,\sigma^2|Y) = \left( \frac{n}{2\pi \mathrm{e}} \right)^{n/2} |(Y-X\gamma)'(I - P)(Y-X\gamma)|^{-n/2} \label{}$

がわかった.

あとは, $P = P_{X|H}$ , $\mathrm{rank} \ P = r - k$ を示せばよい. $P$ と $(G'X'XG)^-(G'X'XG)$ はべき等であること, 及び, 命題4から,

$P$ の定義から明らかに $\mathrm{Im} \ P \subset \mathrm{Im} \ XG$ であることから, $\mathrm{Im}\ P = \mathrm{Im} \ XG$ がわかる. $G$ は $\mathrm{Ker} \ H'$ の基底を並べた行列なので, $\mathrm{Ker} \ H' = \mathrm{Im} \ G$ で, $\mathrm{Im} \ H \subset \mathrm{Im} \ X'$ であるから,

$\mathrm{Ker} \ X = (\mathrm{Im} \ X')^{\perp} \subset (\mathrm{Im} \ H)^{\perp} = \mathrm{Ker} \ H' = \mathrm{Im} \ G. \label{}$

これより, $P = P_{X|H}$ であることが分かる. また, $\mathrm{rank} \ XG = \mathrm{rank} \ G - \dim (\mathrm{Ker} \ X) = (p-k) - (p-r) = r - k$ であるから, $\mathrm{rank} \ P = r -k$ . (証明終わり)

以上の結果から,

$\lambda(Y) = \frac{|(Y-X \gamma)'(I-P_X)(Y - X \gamma)|^{n/2}}{|(Y - X \gamma)'(I-P_{X|H})(Y - X \gamma)|^{n/2}} \label{eq03}$

がわかった. 両辺を $2/n$ 乗し,

$F := \frac{n - r}k\frac{(Y - X \gamma)'(P_X-P_{X|H})(Y - X \gamma)}{(Y - X \gamma)'(I-P_X)(Y - X \gamma)} \label{}$

と定めると,

$(\lambda(Y))^{2/n} = \frac{1}{1 + \frac{k}{n-r}F} \quad \cdots (\ast)$

尤度比検定は $\lambda(Y)$ が小さいときに帰無仮説を棄却する検定であり, これは(*)から, $F$ が大きいときに棄却する検定である. よって, この尤度比検定は $F$ による検定と等価である.

この確率変数 $F$ は前回, 制約条件つき最小二乗値を通じて導いた検定統計量と同じであり, 帰無仮説 $H_0$ の下では自由度 $k$ , $n-r$ の $F$ 分布に従う. よって, 尤度比検定から出発して, 最小二乗値を用いて表される検定統計量 $F$ が自然に導かれることが分かった.

2023-09-08

分散分析(ANOVA)のF検定の導出について１——制約つき最小二乗法を用いたアプローチ

前置き
準備
ガウス・マルコフモデル
線形制約条件つき最小二乗法からのF検定の導出
分散分析のF検定の導出
参考文献

前置き

分散分析(analysis of variance, ANOVA)は複数のグループの平均に統計学的に優位な差があるかを調べる方法である. この文章では分散分析のF検定の検定統計量が数学的にどのように導出されるかを説明する. ここでは, 検定統計量の数学的な導出にフォーカスし, 分散分析自体についての説明はしない. これは, 分散分析を学び始める人にとってはこの内容は過度に詳細だろうし, 一定の理解がある人にとっては一からの説明は恐らく冗長だからだ.

このような古典的な内容をなぜ改めて書くのかを簡単に説明する. 私は最近, 数理統計学を勉強している. 私の読んでいるテキスト[久保川]では¹, 分散分析の概要と、F検定で等平均性をテストする方法が説明されているが, その統計量が何故F分布に従うかは書かれていなかった. そこで, 他の手持ちのテキストやウェブ上で色々と調べてみたが, 検定統計量の導出を解説した文章はあまりなく, あったとしても長大なテキストの一部だったり, 線形回帰モデル(もしくは, ガウス・マルコフモデル)の応用として軽く触れられているだけで詳細は省略されていたりした. 最終的には, [吉田], [ラオ]の中に厳密な導出を見つけることができたが, 前者は洗練されているがやや直感的に理解しにくく, 後者は一部証明が理解できなかった. ともあれ, これらを通じて導出の概要はわかり, [吉田]の方法をやや素朴にした(つまり, 洗練度を減じた)方法で証明を与えることができたので, その内容をまとめておくことにした.

また, ここで説明するアプローチは少し天下り的だが, 尤度比検定を出発点として自然にこのアプローチに至る道筋を以下の記事で説明したので, 興味があれば, 追って読んでほしい.

kihanto.com

今回の内容は, self-contained に説明するのは骨が折れるので, いくつかの命題は証明せずに用いる. といっても, それらはどれも適切な予備知識があれば, 容易に証明は理解できるものばかりである. 具体的には非心カイ二乗分布と分布の特性関数についての基本的な知識があれば十分である. それらを除けば, 必要な道具は基本的には, 学部一回生レベルの線形代数, 微分積分であり, これらは前提知識として断りなく用いる.

ここからの内容のアウトラインを説明する. まず準備として, 線形代数から単純ながらも鍵となる命題を2つ思い出し, 今回の議論で中心的な役割を果たすフィッシャー・コクランの定理を証明なしで述べる. 次に, ガウス・マルコフモデルを導入する. これは分散分析モデルを含む, より一般的なモデルであり, これに対し, 制約つき最小二乗値の従う分布を導出する. 最後に, ガウス・マルコフモデルに対して得た結果を分散分析モデルに適用することで, よく知られているF検定の検定統計量が導出されることを示す.

準備

まず, 線形代数から2つの命題を準備する. いずれもほぼ常識と言える命題である.

1つ目の命題は内積と行列の転置の基本的な関係についてである.

命題1 行列 $A$ に対し, $\mathrm{Ker} \ A' = (\mathrm{Im} \ A)^{\perp}$ .

証明ベクトル $x,y$ に対し, $x \in \mathrm{Ker} \ A'$ なら, $x'Ay = (A'x)'y = 0$ . $y$ は任意に取れるので, $x \in (\mathrm{Im} \ A)^{\perp}$ . 逆に, $x \in (\mathrm{Im} \ A)^{\perp}$ なら, $(A'x)'y = x'Ay = 0$ . $y$ は任意より, $A'x = 0$ . つまり, $x \in \mathrm{Ker} \ A'$ . (証明終わり)

線形変換 $P$ が $P^{2} = P$ , $P' = P$ を満たすとき, 直交射影という. 部分空間 $U$ が $P$ の像であるとき, $P$ を $U$ への射影という. 2つ目の命題は, 直交射影が任意のベクトルを射影の像に最も近い点に移すことを示すものである.

命題2

$P$ を部分空間 $U$ への直交射影とするとき, 任意の $y \in \mathbb{R}^{n}$ に対し, $\min_{z \in U}|y - Pz|^{2} = |y - Py|^{2} = y'(I - P)y.$

証明

任意の $z \in U$ に対し, $|y - Pz|^{2} = |y - Py|^{2} + |P(y-z)|^{2} + 2(y-Py)'P(y-z).$ ここで, $(I-P)P = 0$ より, $(y-Py)'P(y-z) = y'(I-P)P(y-z) = 0$ . よって, $|y - Pz|^{2} \geq |y - Py|^{2}$ . (証明終わり)

次に, フィッシャー・コクランの定理と呼ばれる定理を紹介する. これは分散分析の検定統計量の導出において, 最も基本的な役割を果たす定理で, 単位行列の分解とガウス確率変数の二次形式の関係を示している. 自由度 $n$ , 非心度 $\delta$ の非心カイ二乗分布を $\chi^{2}(n,\delta)$ と書く. 非心カイ二乗分布については例えば, [吉田], Wikipedia 参照. ただし, Wikipedia の日本語版記事はかなり内容が薄いので, 英語版を勧める.

注意: 非心カイ二乗分布はややマイナーであり, また, その確率密度関数には無限和が表れ, 厄介そうに見えるかもしれない. しかし, カイ二乗分布と同様の再生性をもち, 特性関数(もしくはモーメント母関数)もかなり単純で, 意外に扱いやすい対象である.

$Y = (Y_1, \dots, Y_n)'$ を $n$ 次元のガウス確率変数で平均ベクトル $\mu = (\mu_1, \dots,\mu_n)'$ , 分散共分散行列 $\sigma^2 I$ ( $I$ は n次単位行列) をもつものとする.

定理1 (フィッシャー・コクランの定理)

$A_1, \dots, A_k$ をn次対称行列で, $A_1 + \cdots + A_k = I$ を満たすものとする. $Q_i := Y'A_i Y$ , $\mathrm{rank} \ A_i := n_i$ とする. このとき, 以下は同値:

各 $1 \leq i \leq k$ に対し, ある $\delta_i \geq 0$ が存在し, $Q_i \sim \chi^2(n_i,\delta_i)$ かつ $Q_1, \dots, Q_k$ は独立.
$n_1 + \cdots + n_k = n$ .
各 $1 \leq i \leq k$ に対し, $A_i^2 = A_i$ .
各 $1 \leq i,j \leq k$ , $i \neq j$ に対し, $A_i A_j = O$ ( $O$ はゼロ行列).

またこのとき, $\delta_i = \mu' A_i \mu$ , $\mu_1^2 + \cdots + \mu_k^2 = \delta_1 + \cdots + \delta_k$ .

証明は省略する. 例えば, [吉田]参照. 以下で用いるのは, 次の系で与えられる, この定理の特別な場合のみである.

系1

$P_1$ , $P_2$ を $\mathbb{R}^n$ 上の直交射影とする. $Q_i := Y'P_iY \ (i = 1,2)$ , $n_i := \mathrm{rank} \ P_i$ , $\delta_i := \mu'P_i\mu$ とする. このとき, 以下が成立:

$Q_i \sim \chi^2(n_i,\delta_i)$ .
$P_1$ と $P_2$ が互いに直交するとき, つまり, $P_1P_2 = P_2P_1 = O$ のとき, $Q_1$ と $Q_2$ は独立.

証明

(i): 分解 $I = P_i + (I - P_i)$ は明らかに $P_i(I - P_i) = O$ を満たすので, フィッシャー・コクランの定理から $Q_i := Y'P_iY$ がわかる.

(ii): 直交性の仮定から, $(I - P_1 - P_2)^2 = I - P_1 - P_2$ . よって, 分解 $I = P_1 + P_2 + (I-P_1-P_2)$ にフィッシャー・コクランの定理を適用することでわかる. (証明終わり)

この系も含め, 以下では単にフィッシャー・コクランの定理と呼ぶ.

ガウス・マルコフモデル

$Y = X \beta + \varepsilon. \label{}$

ここで,

で, $(\varepsilon_i)_{1 \leq i \leq n}$ は独立同分布で平均 $0$ , 分散 $\sigma^2$ の正規分布に従うとする. $Y$ と $X$ は既知の値として扱い, パラメータ $\beta$ と分散 $\sigma^2$ は未知の値である. $r := \mathrm{rank} \ X \leq n-1$ を仮定する. これは全く制限的な仮定ではない. またいくつかのテキストでは仮定されているようだが, ここでは $X$ がフルランクであること, つまり, $\mathrm{rank}\ X = p$ は仮定しない. これを仮定すると分散分析モデルに適用できないためである.

ガウス・マルコフモデルは分散分析モデルよりも一般的なモデルである. 例えば, 1つの因子が3つの水準をもつ一元配置モデル $Y_{ij} = \alpha_i + \varepsilon_{ij} \quad (1 \leq i \leq 3, \ 1 \leq j \leq 2)$ とするとき,

とすればよい.

$H_0: \ H'\beta = d, \quad H_1: \ H'\beta \neq d \label{}$

線形制約条件 $H'\beta = d$ はパラメータ間に制約を与える. 例えば上の一元配置の例で,

とすると, $H'\beta = d$ は $\alpha_1 = \alpha_2 = \alpha_3$ と同値であり, 3つのパラメータが等しいという制約条件である. よって, $H$ をこのようにとれば, 帰無仮説 $H_0$ , 対立仮説 $H_1$ による検定は一元配置の等平均性の検定と一致する.

線形制約条件つき最小二乗法からのF検定の導出

$H'\beta = d$ を満たす $\beta$ を $1$ つ固定し, $\gamma$ とする. $\beta$ , $\sigma^2$ の真値を $\beta_0$ , $\sigma_0^2$ とする. 次の定理が最も重要である. 制約条件つきの最小二乗値が非心カイ二乗分布に従うことを示している.

定理2

以下を定める:

$R_1^2 := \min_{\beta:H'\beta = d} |Y-X\beta|^2. \label{}$

このとき,

$R_1^2 = (Y-X\gamma)'(I - P_{X|H})(Y-X\gamma). \label{}$

ここで, $P_{X|H}$ は直交射影で, その像は

$\{ X\beta \mid H'\beta = 0 \} \cdots (*)$

で与えられ, $\mathrm{rank} \ P_{X|H} = r-k$ . よって, $R_1^2 / \sigma_0^2$ は非心度 $(X(\beta_0-\gamma))'(I-P_{X|H})(X(\beta_0-\gamma)) / \sigma_0^2$ , 自由度 $n-(r-k)$ の非心カイ二乗分布に従う.

証明

まず, $\mathrm{rank} \ H' = k$ より, $\mathrm{dim} \ (\mathrm{Ker} \ H') = p-k$ である. $H'\gamma = d$ より, $\mathrm{Ker} \ H'$ の基底を $1$ つとり, $g_1,\dots,g_{p-k}$ とし, この基底を並べて作った $p \times (p-k)$ 行列 $G := (g_1, \dots,g_{p-k})$ とすると, $H'\beta = d$ を満たす $\beta$ の集合は $\{ \gamma + G\delta \mid \delta \in \mathbb{R}^{p-k} \}$ である. よって, $R_1^2 = \min_{\delta \in \mathbb{R}^{p-k}}|(Y - X\gamma) - XG\delta|^2$ . 命題2より, $R_1^2 = (Y-X\gamma)'(I - P)(Y - X\gamma)$ . ここで, $P$ は $\mathrm{Im} \ XG$ への直交射影. よって, $\mathrm{Im} \ XG$ が(*)で与えられ, $\mathrm{rank} \ P = r-k$ であることを示せばよい.

$G$ は $\mathrm{Ker} \ H'$ の基底を並べた行列なので, $\mathrm{Ker} \ H' = \mathrm{Im} \ G$ . よって, $\mathrm{Im} \ XG = \mathrm{Im} \ X|_{\mathrm{Ker} \ H'}. \label{}$ ただし, $X|_{\mathrm{Ker} \ H'}$ は線型写像 $X$ の定義域を部分空間 $\mathrm{Ker} \ H'$ に制限したものである. これより, $\mathrm{Im} \ XG$ が(*)に一致することがわかる. $\mathrm{Im} \ H \subset \mathrm{Im} \ X'$ であることを思い出すと,

$\mathrm{Ker} \ X = (\mathrm{Im} \ X')^{\perp} \subset (\mathrm{Im} \ H)^{\perp} = \mathrm{Ker} \ H' = \mathrm{Im} \ G. \label{}$

よって, $\mathrm{rank} \ XG = \mathrm{rank} \ G - \dim (\mathrm{Ker} \ X) = (p-k) - (p-r) = r - k$ . ゆえに, $\mathrm{rank} \ P = r -k$ . 最後の主張はフィッシャー・コクランの定理の直接的な帰結である. (証明終わり)

次に, $R_0^2 := \min_{\beta} |Y-X\beta|^2 \label{}$ と定めると, 命題2より, $R_0^2 = Y'(I-P_X)Y. \label{}$ ただし, $P_X$ は $\mathrm{Im} \ X$ への直交射影. よって, フィッシャー・コクランの定理より, $R_0^2/ \sigma_0^2$ は自由度 $n-r$ のカイ二乗分布に従う. ここで非心度が $0$ となるのは, $(X\beta_0)'(I - P_X)X\beta_0 = 0$ だからである. また明らかに, $R_0^2 = (Y-X\gamma)'(I-P_X)(Y-X\gamma)$ である. これより,

$R_1^2 - R_0^2 = (Y-X\gamma)'(P_X - P_{X|H})(Y-X\gamma) \label{}$

である. 定義から明らかに, $\mathrm{Im} \ P_{X|H} \subset \mathrm{Im} \ P_X$ であるから, $(R_1^2 - R_0^2) / \sigma_0^2$ は自由度 $r - (r-k) = k$ , 非心度 $(X(\beta_0-\gamma))'(P_X-P_{X|H})(X(\beta_0-\gamma)) / \sigma_0^2$ の非心カイ二乗分布に従う. また, $P_X - P_{X|H}$ と $P_{X|H}$ は互いに直交する直交射影であるから, フィッシャー・コクランの定理より, $R_1^2 - R_0^2$ と $R_0^2$ は独立である. 以上より,

$F := \frac{n-r}k\frac{R_1^2 - R_0^2}{R_0^2} \label{}$

は自由度 $k$ , $n-r$ , 非心度 $(X(\beta_0-\gamma))'(P_X-P_{X|H})(X(\beta_0-\gamma)) / \sigma_0^2$ の非心 $F$ 分布に従うことがわかった. 非心F分布については例えば, 英語版のWikipedia 参照 (日本語版記事は無いようである).

これで帰無仮説 $H_0$ , 対立仮説 $H_1$ の場合の検定統計量 $F$ が構成できたように見えるが, $F$ の非心度には知ることのできない分散の真値 $\sigma_0^2$ が含まれている. しかし, これは帰無仮説の下では $0$ になるので問題はない.

命題3

$\beta$ の真値 $\beta_0$ が帰無仮説 $H_0$ を満たすとする. すなわち, $H'\beta_0 = d$ とする. このとき,

$(X(\beta_0-\gamma))'(P_X-P_{X|H})(X(\beta_0-\gamma)) = 0. \label{}$

証明

$H'\beta_0 = H'\gamma = d$ より, $\beta_0 - \gamma \in \mathrm{Ker} \ H'$ . よって, $X(\beta_0-\gamma) \in \mathrm{Im} \ P_{X|H}$ . (証明終わり)

よって, 帰無仮説の下では $F$ は非心度 $0$ の通常の $F$ 分布に従うことがわかり, $F$ は $X$ と $Y$ のみから計算できるので, これで帰無仮説 $H_0$ , 対立仮説 $H_1$ の検定を行うことができる.

分散分析のF検定の導出

前節までの結果を用いて, 二元配置の分散分析モデルで交互作用のない場合の検定量を導出する. より複雑な分散分析モデルでも, 計算が多少複雑になるだけで同様に導出できる. モデルは以下で与えられる:

$Y_{ij} = \alpha_i + \beta_i + \varepsilon_{ij} \quad (1 \leq i \leq I, 1 \leq j \leq J). \label{}$

これがどのようなガウス・マルコフモデルに対応するかを念のため述べておく. 例えば, $I = 2$ , $J = 3$ のときは,

と定めればよい. このように定めるとき, 掃き出し法により, $\mathrm{rank} \ X = I+J-1$ がわかる(やや計算と考察を要する).²

このとき, 以下の帰無仮説 $H_0$ , 対立仮説 $H_1$ の検定を考える:

$H_0: \ \alpha_1 = \cdots = \alpha_I =:\alpha, \quad H_1: \ \alpha_i \neq \alpha_k \quad (\exists \ i,k). \label{}$

既にみたように帰無仮説 $H_0$ は適当な行列 $H$ を用いた線形制約条件で表現することができる. また, 容易に $\mathrm{rank} \ H = I-1$ がわかる.

前節の結果より, 検定を行うには,

$R_0^2 = \min_{\alpha_1,\dots,\alpha_I,\beta_1,\dots,\beta_J}\sum_{i,j}(Y_{ij} - \alpha_i - \beta_j)^2, \quad R_1^2 = \min_{\alpha,\beta_1,\dots,\beta_J}\sum_{i,j}(Y_{ij} - \alpha - \beta_j)^2 \label{}$

を求め, 自由度 $I-1$ , $IJ - I - J - 1 = (I-1)(J-1)$ の $F$ 分布に従う確率変数

$F = \frac{(I-1)(J-1)}{I-1}\frac{R_0^2 - R_1^2}{R_0^2}$

の値が棄却域に入るか否かで検定を行えばよい. $R_0^2$ , $R_1^2$ の値は, 例えば偏微分するなどの通常の方法で求めることができ,

$R_0^2 = \sum_{i,j}(Y_{ij} - \bar{Y}_{\cdot j} - \bar{Y}_{i \cdot} + \bar{Y}_{\cdot \cdot})^2, \quad R_1^2 = \sum_{i,j}(Y_{ij} - \bar{Y}_{\cdot j})^2 \label{}$

がわかる. ただしここで,

$\bar{Y}_{\cdot j} := \frac{\sum_iY_{ij}}I, \quad \bar{Y}_{i \cdot} := \frac{\sum_jY_{ij}}J, \quad \bar{Y}_{\cdot \cdot} := \frac{\sum_{i,j}Y_{ij}}{IJ}.$

また容易に以下がわかる:

$R_0^2 = \sum_{i,j}( (Y_{ij} - \bar{Y}_{\cdot j}) - (\bar{Y}_{i \cdot} - \bar{Y}_{\cdot \cdot} ) )^2 = R_1^2 - J\sum_i(\bar{Y}_{i \cdot} - \bar{Y}_{\cdot \cdot})^2. \label{}$

よって,

$R_1^2 - R_0^2 = J \sum_i(\bar{Y}_{i \cdot} - \bar{Y}_{\cdot \cdot})^2. \label{}$

以上より,

$F = \frac{(I-1)(J-1)}{I-1}\frac{J\sum_i(\bar{Y}_{i \cdot} - \bar{Y}_{\cdot \cdot})^2}{\sum_{i,j}(Y_{ij} - \bar{Y}_{\cdot j} - \bar{Y}_{i \cdot} + \bar{Y}_{\cdot \cdot})^2}$

これは分散分析において通常用いられている検定量であるから, これで分散分析の検定統計量が $F$ 分布に従うことが証明できた.

最後にややどうでもいいコメントをしておくと, 交互作用のない二元配置ぐらいの比較的単純なモデルならば, $R_0^2$ と $R_1^2$ の式を展開し, $R_0^2 = Y'PY$ , $R_1^2 = Y'QY$ なる直交射影 $P,Q$ を具体的に求め, それらに対してフィッシャー・コクランの定理を適用することで, 確率変数 $F$ が $F$ 分布に従うことを示すことができる. ただし, このような力技はモデルが複雑になるにつれて破綻するように思える.³ 一方で, 今回の方法はモデルが複雑化してもそこまで難しくはならないし, 検定統計量の具体系が予めわかっていなくても, 最小値問題を解くことで機械的に検定量を導出できるという点で利点がある.

参考文献

・[久保川]: 久保川達也. 現代数理統計学の基礎. 共立出版, 2017.

・[吉田]: 吉田朋広. 数理統計学. 朝倉書店, 2006.

・[ラオ]: ラオ, Ｃ. ラダクリシュナ. 統計的推測とその応用. 東京図書, 1986.

参考文献リストは文章末尾にある。↩
このランクの計算がより一般のモデルの場合に容易にできるかは分からない.↩
ただしある種の組合せ論的な手法と合わせれば, 一般の場合に通用する議論ができそうな気はしている.↩

2023-08-24

ロト6の数字の偏りについて

1. はじめに
2. 目的
3. 方法
4. 結果
5. 結論
6. 余談
参考文献
付録: 検定に用いたコード

1. はじめに

この文章ではロト6の当選数字の出方に偏りがあるかどうかを統計的に考える。

ロト6とは、多くの人は知っていると思うが宝くじの一種で、1から43までの数字から6個の数字を選び、6個の当選番号と3個以上一致すると、一致した数字の個数に応じて賞金がもらえるというものである。詳しいルールについては以下のオフィシャルな説明を読んでほしい。

www.takarakuji-official.jp

考えたいのは、この数字の出方に統計的な偏りがあるかである。ここでは「統計的」という部分に力点があり、単に偏りがあるかという話ではない。

例えば、ロト6の第1回から第1820回までで最もよく出た数字は6で282回、最も出ていない数字は9で225回であり、57回分の差がある。しかしこれをもって、ロト6では6を選ぶと当たりやすく、9を選ぶと外しやすいと結論することはできない。全ての数字が同じ確率で現れるとしても、偏りは生じるからだ。

偏りには自然なものと不自然なものがある。例えば、コインを100回投げて表が20回しか出なかったら、裏が出やすいコインだと疑うのは当然である。しかし、表が35回出るのが自然か不自然かというと、判断に困るのではないだろうか。同じく、ロト6の57回分の差も、これが自然かどうかは直感的には判断しにくい。

統計学を用いると、自然な偏りと不自然な偏りを区別できる。特に適合度検定(の一般化)を使って、ロト6の数字の出方に偏りがあるかを考えるのが、この文章の趣旨である。ただし、統計学の結論は確率に基づいているため100%確かではなく、常に一定の不確実性があることは理解しておいてほしい。

2. 目的

過去のロト6の結果(第1回から第1820回)までを集計し、各数字および2つの数字の組の出現頻度に偏りがあるかを調べる。これは例えば、特定の数字(例えば5)や数字の組(例えば11と13の組)が当選番号に出やすいというような偏りがあるかを調べるということ。また、3つ以上の数字の組の出やすさは調べない。これは今回の手法では調べられないからである(以下の「方法」のセクションの末尾の注意を参照)。

注意: 頻度のみに着目して調べるため、抽選同士の相関は分からない。例えば、ある回とその次の回の当選番号に相関があるかは今回の考察からはわからない。

3. 方法

(注意: 統計学をあまり知らなかったり、数式に抵抗があったりする人にとってはこれ以降は読みづらいと思うが、細かい部分は飛ばしても大まかな論旨は把握できるように書いたつもりなので、適当に飛ばし飛ばし読んでほしい。)

よく知られた適合度検定を修正したものを用いる。適合度検定は統計学を学んだ人なら誰でも知っている方法だが、今回のケースには素朴には適用できず、修正が必要になる。

3.1. 適合度検定とは

念のため、通常の適合度検定について簡単に振り返る。詳しくは以下のWikipediaの記事または統計学の入門書を参照。

ja.wikipedia.org

ある試行において、 $E_1, E_2, \dots, E_N$ の $N$ 個の結果が起こりうるとする。これらの結果は1回の試行でただ1つだけが起こるとする。それぞれの結果が起こる確率を $p_1, p_2, \dots, p_N \ (p_i \geq 0, p_1 + p_2 + \cdots + p_N = 1)$ と予想したとき、これが正しいか否かを検定するのが適合度検定である。¹

$n$ 個の観測結果 $Y_1, Y_2, \dots, Y_n$ が与えられたとしよう。つまり、 $Y_i$ はそれぞれ $E_1, E_2, \dots, E_n$ のいずれかである。このとき、検定統計量は

$X^{2} = \sum_{i=1}^{n}\frac{(O_{i}-np_{i})^2}{np_i}$

である。ここで、 $O_i$ は結果 $E_i$ の観測度数、つまり、 $Y_i = E_i$ となるような $i$ の個数である。 $X^{2}$ の値は $O_i$ が期待される観測回数 $np_i$ から離れるほど大きくなる。よって、 $X^{2}$ が大きいほど予想が外れている可能性が高いことになる。

$X^{2}$ の分布は $n \to \infty$ の極限で、パラメータ $N-1$ のカイ二乗分布に近づくため、大きい $n$ について $X^{2}$ がカイ二乗分布に従うと見なすことでカイ二乗分布表を用いて検定を行うことができる。 $n$ をどれくらい大きくとればいいかについては実用的な基準として、すべての $i$ に対し $np_i \geq 10$ となれば良いとされているようである。

適合度検定のイメージを掴んでもらうため、簡単な例を1つ挙げる。サイコロの各目が同じ確率で出るかを検定することを考える。各目 $i \ (1 \leq i \leq 6)$ が出るという結果を $E_i$ とし、 $p_i=1/6$ と定める。サイコロを $n$ 回振って、上に挙げたカイ二乗値 $X^2$ を計算する。この値は各目が出る頻度が $1/6$ からズレるほど大きくなる。よって、これがある値よりも大きい場合、 $p_i=1/6$ という仮説は棄却され、サイコロは公正ではないと見なされる。

3.2. 適合度検定の修正の必要性

ロトで各数字の出現確率を検定するには、上に書いたような適合度検定は適用できず、工夫が必要となる。

素朴に考えると、1回の抽選で6個選ばれた数字のうちに $i$ が含まれることを結果 $E_i$ として適合度検定を行えばいいと思うかもしれない。しかし、1回の抽選で複数の数字が同時に選ばれるため、これでは異なる $i$ と $j$ に対して $E_i$ と $E_j$ は同時に起こらないという適合度検定の前提に反してしまう。

では、6個の数字を選ぶことを、1個の数字を6回選ぶことと見なせばいいと思うかもしれないが、これにも問題がある。なぜなら、各数字の選択確率は、1回目から6回目まで変わってしまうからだ。例えば、1回目に取り出した数字は次の5回では選ばれない。このことは、観測結果が独立でなければならないという適合度検定の前提(上では明示的には書かなかったが)に反してしまう。

結局、ロトに関しては、単純な適合度検定を適用することはできない。検定統計量の導出から修正を加える必要がある。幸いそのような場合に適合度検定を一般化した先行研究がある: Joe (1993), Genest et. al. (2002). ²

これら2つの先行研究ではそれぞれ、ロトのように複数の数字からいくつかの数字を同時に選出する試行に対し、数字(または数字の組)の出現確率が等しいかを検定する方法を与えている。この2つの論文では異なる形の検定統計量を与えているが、各数字(つまり1つの数字)の出現確率の場合は統計量は等しくなる。

今回は特に理由はないが、Joe (1993) ではなく Genest et. al. (2002) の方法を採用する。

3.3. 検定統計量

以下では観測回数を $n$ と表している。今回の場合は1820回分の結果がわかっているので $n=1820$ である。 Genest et. al. (2002)の結果から、以下のように検定統計量を定めることができる。

各数字の出現確率の等しさの検定

検定統計量を次のように定める:

$U_1 = \frac{42}{37}\sum_{i=1}^{43} \frac{(O_i - 6n/43)^{2}}{(6n/43)}$ .

$U_1$ の分布は $n \to \infty$ の極限において、パラメータ42のカイ二乗分布に近づく。

$\frac{6}{43} = \binom{42}{5} / \binom{43}{6}$ は各数字の理論出現確率 (選出に偏りがない場合に、1回の抽選で特定の数が出る確率)。
よく見ると、 $U_1$ は上で言及した、数字を1つずつ選出したとみなして計算した場合のカイ二乗値の42/37倍。よって、通常の適合度検定の計算を流用できる。

2つの数字の組の出現確率の等しさの検定

検定統計量を次のように定める:

$U_2 = \sum_{i = 1}^{43}\sum_{j = i+1}^{43}\frac{(O_{i,j} - E_{i,j})^{2}}{E_{i,j}}$ .

$E_{i,j} := n \frac{5 \cdot 6}{42 \cdot 43}$ で、これは各ペアの理論出現確率のn倍。
$O_{i,j}$ はペア $(i,j)$ の観測度数。

$U_2$ の分布は $n \to \infty$ の極限において、次の確率変数の分布に近づく: $w_1 K_{42} + w_2 K_{860}$ .

ここで、 $K_{42}$ と $K_{860}$ はそれぞれパラメータ $42$ と $860$ のカイ二乗分布に従う独立な確率変数であり、 $w_1 = (5 \cdot37) / 41$ 、 $w_2 = (36 \cdot 37) / (40 \cdot 41)$ 。

上記2つの統計量の極限分布について、数値計算によってp値を求めた。

注意: 3つ以上の数字の組について検定を行わない理由

先に言及した2つの先行研究では3つ以上の数字の組に対しても検定統計量を与えている。しかし、それらに対して検定を行うには現状ではデータが少なすぎると考えられる。既に述べたように、適合度検定を適用する実用的な基準として、(観測の総回数) × (理論出現確率) が10以上というのがある。この基準が今回の修正した適合度検定に対しても適切かは明らかではないが、ここでは同じ基準を適用することにする。このとき、それぞれの場合について(観測の総回数) × (理論出現確率)の値を求めると以下のようになる:

1つの数字のときは、 $1820 \cdot (6/43) \fallingdotseq 254$ .

2つの数字の組のときは、 $1820 \cdot (5 \cdot 6 ) / (42 \cdot 43) \fallingdotseq 30$ .

3つの数字の組のときは、 $1820 \cdot (4 \cdot 5 \cdot 6 ) / (41 \cdot 42 \cdot 43) \fallingdotseq 3$ .

このように3つ以上の数字の組に対しては10を下回ってしまう。よって、現状ではこの検定統計量を用いるのは恐らく適切ではない。

4. 結果

以下が結果である。

各数字の出現確率の等しさの検定のp値 ≒ 0.88.
2つの数字の組の出現確率の等しさの検定のp値 ≒ 0.87.

5. 結論

初めに少しだけp値についての一般的な説明をする。 p値は0から1の間の値を取り、この値が小さいとき、仮説 (今回の場合は例えば全ての数や組が等確率で現れるという仮説)は怪しいとされる。通常、0.05や0.01などのしきい値を設定し、これより小さい場合は統計的に有意な差があるとされ、元々の仮説は棄却される。注意点としては、p値が大きいほど仮説の信頼性が一層高まるわけではなく、単にデータが仮説と整合していることを示しているに過ぎない事である。

結論を言うと、今回のp値は極めて大きい。よって、等確率性を疑うことは難しい。

6. 余談

先行研究である Genest et. al. (2002) ではカナダのロトの1798回分の結果に対し検定を行い、p値はそれぞれ0.104, 0.300としている。これはそれなりに大きいとはいえ、日本のものほど極端に大きくはない。カナダのロトは49個の数字から6個を選ぶので、日本のものとは少し異なるが、それだけでこれほど差が出るとは思えない。また、同論文ではボーナスナンバーを含めた場合、各数字の出現確率の検定のp値は0.044となると報告している。これは有意水準5%の場合には仮説は棄却されてしまうほど低い。カナダのロトの等確率性はやや怪しいのかもしれない。³

ついでなので、日本のロトの場合も、ボーナスナンバーを含めた等確率性の検定を各数字、2つの数字の組にそれぞれ試してみたところ、どちらもp値はほぼ0.90となった。やはりp値は十分に大きい。

参考文献

・H. Joe. Tests of uniformity for sets of lotto numbers. Statist. Probab. Lett., 16(3):181–188, 1993.

・C. Genest, R. A. Lockhart, and M. A. Stephens. $\chi^{2}$ and the lottery. The Statistician, 51(2):243–257, 2002.

付録: 検定に用いたコード

検定に用いたコードを載せておく。私は簡単なコードを読めるぐらいの実力しかないので、 ChatGPTに大いに頼ってプログラムを組んだ。多分正しく動作していると思うが確証はない。ロトの当選番号のデータ(コード内では"loto_results.txt"として読み込まれている)はここには載せないが簡単に手に入る。

各数字の出現確率の等しさの検定

コードを開く

# 各数字の等確率性の検定
import numpy as np
from scipy.stats import chisquare, chi2

# 1からlまでの数値
l = 43

with open('loto_results.txt') as f:
    nums = np.array([int(x) for x in f.read().split()])

# ビンのカウント
bins, _ = np.histogram(nums, bins=range(1, l+2))

# 1からlまでの数値が等確率で選ばれるという仮説に基づいて、期待される頻度を計算
expected_frequencies = [len(nums) / l] * l

# カイ二乗適合度検定を実行
chi2_stat, p_val = chisquare(bins, expected_frequencies)

# カイ二乗統計量の調整
adjusted_chi2_stat = chi2_stat * (42 / 37)

# 自由度42のカイ二乗分布でp値を計算
adjusted_p_val = 1 - chi2.cdf(adjusted_chi2_stat, 42)

print("統計量:", adjusted_chi2_stat)
print("p値:", adjusted_p_val)

統計量: 31.529937629937628 p値: 0.8809230858481042

2つの数字の組の出現確率の等しさの検定

統計量の数値積分がうまくできなかったので、モンテカルロ法を用いて積分の近似値を求め、そこからp値を算出している。

コードを開く

# 2つの数の組の等確率性の検定
import numpy as np
from itertools import combinations
from scipy.stats import chi2

# 43 x 43の全ての成分が0からなる配列pairsを作る
pairs = np.zeros((43, 43), dtype=int)

# ファイルから読み取り
with open('loto_results.txt') as f:
    given_array = np.array([int(x) for x in f.read().split()])

# 6個ずつに区切る
nums = np.array_split(given_array, len(given_array) // 6)

# numsの各要素xの2つの要素の組(i,j)に対し、pairs[i][j] += 1する
for num_group in nums:
    for i, j in combinations(num_group, 2):
        small, large = min(i, j), max(i, j)
        pairs[small-1][large-1] += 1

# 定数E
E = 1820 * 5 * 6 / (42 * 43)

# (pair[i][j] - E)^{2} / E をi < jで和を取る計算
result = sum((pairs[i][j] - E)**2 / E for i in range(pairs.shape[0]) for j in range(i+1, pairs.shape[1]))

# モンテカルロ法によるp値の計算

# 定数とパラメータ
a = 5 * 37 / 41  
b = 36 * 37 / (40 * 41)
df_U = 42
df_V = 860
p = 0.1
n_samples = 10000000  # サンプル数

# UとVのサンプルを生成
samples_U = chi2.rvs(df_U, size=n_samples)
samples_V = chi2.rvs(df_V, size=n_samples)

# Wのサンプルを計算
samples_W = a * samples_U + b * samples_V

# サンプル集合samples_Wと値xが与えられたと仮定
p_value = np.mean(samples_W > result)

print("統計量:", result)
print("p値:", p_value)

統計量: 828.0200000000001 p値: 0.8728807

「正しいか否か」という言い方は厳密には正しくないが、あまり細かいことは気にしないことにする。↩
これらの論文の計算の細部は(少なくとも私にとっては)巧妙であり、完全にはフォローしていない。ただ証明のアイディア自体は通常の適合度検定と同じである。つまり、各数字(または数字の組)の出現回数の共分散行列を求め、それをスペクトル分解することである。↩
カナダ、日本のロトでは実は6個の当選番号だけでなく、追加で1個のボーナスナンバーが選ばれ、ボーナスナンバーはやや変則的な形で当選に関わる。詳細は上記のオフィシャルサイトを参照。↩

2023-08-15

「ピータースのコイントス」について

はじめに
ピータースのコイントス
パラドキシカルな点
メカニズム：大雑把な説明
メカニズム：正確な説明

はじめに

この文章では以下の記事(以下では「元記事」と呼ぶ)について解説する。

gigazine.net

この解説では、元記事は読んでいるものとして進めるが、念のため次のセクションで簡単に要点を振り返る。

元記事では、ピータースのコイントスという賭けのモデルについて、シミュレーションの結果をいくつか並べ、何かしらパラドキシカルなことが起こっていることを示している。しかし、それが結局どういうメカニズムで生じているかは説明していない。この文章ではその部分を整理し説明する。メカニズムを数学的に厳密に説明するには大学レベルの確率論の知識が必要になる。恐らくこれは多くの人にはハードルが高いため、ここでは、

それほど前提知識を必要としない大雑把でやや非厳密な説明
数学的に厳密な説明

の2つに分けて説明する。

ピータースのコイントス

次の賭けを考える。

最初は1ドルを持った状態からスタートする。¹
公平なコイン、つまり表裏がそれぞれ1/2で出るコインを投げ、表が出れば、所持金が1.5倍に、裏が出れば所持金は0.6倍になる。
賭けは好きな回数行うことができる。

問題とするのは、この賭けに参加すべきか、すべきでないかである。そしてこれを素朴に考えると、パラドキシカルなことが起こるのを示していたのが元記事である。その内容を振り返る。

パラドキシカルな点

賭けへの参加・不参加の判断には何らかの基準が必要で、その基準として「期待値」は最も一般的に知られたものだろう。そこで期待値を計算する。確率変数 $X$ は $\mathbb{P}(X = 1.5) = \mathbb{P}(X = 0.6) = 1/2$ を満たすとする。これの期待値は $\mathbb{E}(X) = \frac{1}{2} \cdot 1.5 + \frac{1}{2} \cdot 0.6 = 1.05$ になる。
これは1を超えているから、期待値から見ると有利な賭で、1回の賭けで平均的には所持金は1.05倍になることが期待できそうである。

しかし、元記事の中のシミュレーションで示されていたように、賭けを何回も繰り返すと所持金は減っていく。この記事の後半で説明するように、実は賭けを繰り返すと確実に所持金はどこまでも0に近づいてしまう。 ここにピータースのコイントスの奇妙な点がある。つまり、期待値でみると所持金は増大しそうなのに、実際には所持金は減ってしまう。

ただし一方で、所持金が減っていくのは当たり前という見方もある。それは、この賭けでは1度勝って、1度負けると、所持金は $1.5 \cdot 0.6 = 0.9$ 倍になるからである。この賭けは公平なコインで決まるから、何度も繰り返せば勝率は大体1/2に近づく。だから大まかにいって、1回賭けをする毎に所持金は $\sqrt{0.9}$ 倍になり、所持金が減っていくことはこれで説明がつく。以上から、次のような結論を導けそうである：

この賭けは期待値でみると、有利な賭けに見えるけれども実は不利な賭けである。また、この賭けに対して期待値で有利・不利を判断するのは間違っている。

しかし、こう結論するのは早計である。というのは元記事の別のシミュレーションでは、賭けの参加者数を増やして賭けを何度も行うと、参加者の平均所持金は増えていくことが示されている。この記事の後半で説明するように、実は平均所持金は各回の賭けで概ね $1.05$ 倍となる。この1.05という数字はもちろん先ほど計算した1回の賭けの期待値に由来するものである。

よって、賭けの参加者数を増やした場合には、期待値を基準として賭けが有利だと判断することは正しい。

これらの結果に基づき、元記事では、期待値を基準とした判断は、多数の参加者数の平均を考えるには適切だが、個々人に対しては不適切な場合があることが説明されていた。

なぜ、このようなことが起こるのか？

以上で説明した内容から、以下の2つの観察が得られている：

1人で賭けを繰り返すと、所持金は確実に0に近づく --- (観察A)
参加者が多い場合、平均所持金は概ね $1.05$ 倍のスピードで増えていく --- (観察B)

しかし、この2つは両立しえない。なぜなら、賭けを繰り返すとどの参加者も所持金が0に近づくのなら、どれだけ参加者数を増やしても、その全員の所持金が0に近づき、結局、平均所持金も0に近づくはずだからである。

というわけで、少なくともどちらかの観察は間違っている。ただしこれら2つの観察は、一定の近似のもとでは両立し、元記事のシミュレーションはそのような場合を扱っているため、この相反する事象が観察された。その近似とは何か、そしてその近似の下では何が成り立っているかを説明する。

メカニズム：大雑把な説明

ここではまず大雑把な説明をする。厳密に説明するとややこしくなる部分は少し誤魔化している。まず観察AとBは両立しえないのだから、どちらかが間違っている。実は、観察Bが間違っている。

ただし完全に間違っているわけでなく、次のように修正すれば正しい。

参加者が多い場合、賭けの回数がある程度小さい間は、平均所持金が概ね $1.05$ 倍のスピードで増える確率が高い --- (観察B')

また、(観察A)は正しい。理由は既に上に述べたように、賭けを繰り返すと、大体各回毎に所持金は $\sqrt{0.9}$ 倍になるからである。数学的に厳密な説明は次のセクションで行う。

(観察B')が正しい理由は、参加者数を増やすことで、賭けの後の平均所持金が大数の法則によって期待値に近づくからである。もう少し具体的に言うと、参加者がN人いるとして、そのN人のn回目の賭けの後の所持金をそれぞれ $M_n^{(1)}, M_n^{(2)}, \dots, M_n^{(N)}$ とする。このとき、n回目の賭けの後の平均所持金は $\frac{ M_n^{(1)} + M_n^{(2)} + \cdots + M_n^{(N)} }N$ である。大数の法則より、参加者数Nを無限大とした極限をとると、n回目の賭けの後の所持金の期待値に収束する。つまり、次が成り立つ：

$\lim_{N \to \infty} \frac{ M_n^{(1)} + M_n^{(2)} + \cdots + M_n^{(N)} }N = (1.05)^n.$

この式が示しているのは、参加者数Nを大きくすればn回目の賭けの後の平均所持金は大体 $(1.05)^n$ となるということ。これは(観察B')が成り立つことを示している。

なぜ「賭けの回数がある程度小さい間」という制限が必要かというと、大数の法則による近似を使うためにはNが大きくなければならず、どれだけNを大きくとる必要があるかはnによって異なるからである。 nが大きくなれば、 $M_n^{(i)}$ の分散も大きくなるため、Nもそれに応じて大きくしなくてはならない。そのため、参加者数を事前に大きくとっておいても、ある回数以降の賭けでは近似は機能しなくなる。よって、賭けの回数に制限がつく。

メカニズム：正確な説明

ここでは、上で説明した内容を数学的に説明していく。大学レベルの確率論の知識が必要となる。

まず、賭けの各回の結果を表す確率変数列を定める。 $X_1,X_2, …$ を独立同分布で $P(X_1 = 1.5) = P(X_1 = 0.6) = 1/2$ となる確率変数列とする。
そして、 $M_n = X_1 X_2 \cdots X_n$ と定めると、 $M_n$ は目的通り、賭けのn回目の後の所持金を表す確率変数となることがわかる。

まず観察Aを証明する。 $Y_n := log X_n$ と定めると、やはり $Y_1, Y_2, …$ は独立同分布な確率変数列であり、 $S_n := log M_n$ と定めると、 $S_n = Y_1 + Y_2 + \cdots + Y_n$ となる。つまり、 $S_n$ は独立同分布な確率変数のn個の和なので、ランダムウォークである。1ステップの期待値は $E Y_1 = (1/2) log 1.5 + (1/2) log 0.6 = (1/2) log 0.9$ であり、これは負の数である。よって、大数の強法則から確率1で、 $\lim_{n \to \infty} S_n / n = (1/2) \log 0.9$ 。特に、

$\lim_{n → ∞} S_n = -∞$ .

ここで、 $M_n = \mathrm{e}^{ S_n}$ であることを思い出すと、これは確率1で、 $lim_{ n \to \infty} M_n = 0$ となることを示している。よって、観察Aが示せた。

次に観察B'を示す。参加者が $N$ 人いるとして、各nに対し、 $M_n$ の独立なN個のコピー $M_n^{(1)}, \dots, M_n^{(N)}$ をとる。 n回目の賭けの後の平均所持金 $\bar{M}_{n} := \frac{M_n^{(1)} + \cdots + M_n^{(N)}}N$ と定める。 $\varepsilon, \delta > 0$ を小さくとる。大数の弱法則より、各nに対し、 $N$ によらないある自然数 $k_n$ が存在し、 $N > k_n$ なら、

$\mathbb{P}\left( | \bar{M}_n - (1.05)^n| \leq \varepsilon \right) \geq 1 - \delta.$

これが示していることを解釈すると、以下のようになる：

参加者数を多くすると(Nが大きく)、賭けの回数が小さい間は( $N > k_n$ ならば)、高い確率で(確率 $1-\delta$ 以上で)、平均所持金は期待値に近い(平均所持金と期待値との差は $\varepsilon$ 以下である)。これは観察B'に他ならない。

元記事では100ドルとなっているが、表記の単純化のため変更した。内容には全く影響しない。↩

帰帆島

コメント、ディスカッション歓迎

自転車の走行中の適切な後方確認頻度について

前置き

適切な後方確認頻度とは

具体例における計算

実践と反省

付録: コード

分散分析(ANOVA)のF検定の導出について２——尤度比検定との等価性

前置き

準備

ガウス・マルコフモデル

尤度比検定とF検定の等価性

分散分析(ANOVA)のF検定の導出について１——制約つき最小二乗法を用いたアプローチ

前置き

準備

ガウス・マルコフモデル

線形制約条件つき最小二乗法からのF検定の導出

分散分析のF検定の導出

参考文献

ロト6の数字の偏りについて

1. はじめに

2. 目的

3. 方法

3.1. 適合度検定とは

3.2. 適合度検定の修正の必要性

3.3. 検定統計量

4. 結果

5. 結論

6. 余談

参考文献

付録: 検定に用いたコード

「ピータースのコイントス」について

はじめに

ピータースのコイントス

パラドキシカルな点

メカニズム：大雑把な説明

メカニズム：正確な説明