Weblog on mebius.tokaichiba.jp

中心極限定理の証明の理解でつまずいた

2年前に統計検定準1級の勉強をしていた時、どこかで特性関数を用いた中心極限定理の証明を読んで、なるほど、と思ってノートに書き残した。
先週、そのノートを読んだら、理解できない部分があって、結構な時間悩んでしまった。
同じように悩む人はあまり居ないと思うが、自分用にメモしておく。

中心極限定理とは、期待値がμ、分散がσ²である任意の独立同一分布(i.i.d)に従うn個の確率変数X_iの平均が、その分布がどういう形状であっても、nを大きくすると近似的に期待値μ、分散σ²/nの正規分布に従うというものである。

中心極限定理
$\lim_{n \to \infty} P\left(\frac{\overline{X}-\mu}{\sigma / \sqrt{n}} \le x \right) = \int_{-\infty}^{x} \frac{1}{\sqrt{2\pi}} e^{-\frac{t^2}{2}} dt$

証明方法としては、特性関数の連続定理により、ある確率変数X_nの特性関数φ_n(t)がn→∞とすると全てのtでXの特性関数φ(t)に一致するなら、X_nはXに分布収束するというのを使う。
つまり、X_iの平均を平均μ、分散σ²/nで標準化したZ_nの特性関数φ_Zn(t)が、標準正規分布の特性関数φ_Z(t)に一致することを示せば良い。

標準正規分布に従うZの特性関数は
$\begin{align} \phi_Z(t) & = E[e^{itZ}] = \int_{-\infty}^{\infty} \frac{1}{\sqrt{2\pi}} e^{-\frac{z^2}{2}} e^{itz} dz \\ &= \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty} e^{-\frac{1}{2}((z-it)^2 + t^2)} dz \\ &= \frac{1}{\sqrt{2\pi}} e^{-\frac{t^2}{2}} \int_{-\infty}^{\infty} e^{-\frac{(z-it)^2}{2}} dz \\ &= e^{-\frac{t^2}{2}} \end{align}$
である。最後は
$\int_{-\infty}^{\infty} e^{-a x^2} dx = \sqrt{\frac{\pi}{a}}$
であることを用いた。

本題の、Znの特性関数を計算する。
$\begin{align} \phi_{Z_n}(t) & = E[e^{itZ_n}] = E\left[\exp \left( it\frac{\sqrt{n}}{\sigma}\frac{1}{n}\sum_{j=1}^n (X_j-\mu) \right)\right] \\ &= E\left[\prod_{j=1}^n \exp\left(it\frac{(X_j-\mu)}{\sqrt{n} \sigma} \right)\right] \\ &= \left( E\left[ \exp\left(it\frac{(X_1-\mu)}{\sqrt{n} \sigma} \right)\right] \right)^n \\ \end{align}$

今回、この最後の変形が何故成り立つのかがわからず、調べまくってしまった。結局、九州大学の原先生の公開講座資料の第26頁の「最後のところでは，Xi の分布が同分布であることを用いました」を見て、X_iがi.i.dであることを思い出して解決した。
上の数式では敢えてX₁と書いた（原先生の資料でもそのようになっている）が、筆者のノートでは（おそらくその元資料でも）この部分がXと書かれていたのが、筆者にとって混乱の元だった。

計算を続ける。マクローリン展開により、
$e^x = 1+x+\frac{x^2}{2!}+\frac{x^3}{3!}+\cdots$
なので、2乗の項までを展開して、
$\begin{align} \phi_{Z_n}(t) & = \left( E\left[1 + \frac{it}{\sqrt{n}}\left(\frac{X_1-\mu}{ \sigma}\right) - \frac{t^2}{2n}\left(\frac{X_1-\mu}{ \sigma}\right)^2 - O(\frac{1}{n\sqrt{n}}) \right]\right)^n \\ &= \left(1 - \frac{t^2}{2n} - O(\frac{1}{n\sqrt{n}}) \right)^n \end{align}$
と変形できる。tの1乗の項はE[X_i-μ]=0、2乗の項はE[(X_i-μ)²]=σ²であることを用いた。O(n)は「ランダウの記号」というもので、nを大きくした時に高々f(n)の定数倍となることを表す。3乗の項以降はnを大きくした時に0に収束するので、このように省略している。

$\lim_{n\to\infty} \left(1 + \frac{x}{n} \right)^n = e^x$
なので、
$\begin{align} \lim_{n\to\infty} \phi_{Z_n}(t) & = \lim_{n\to\infty} \left(1 - \frac{t^2}{2n} - O(\frac{1}{n\sqrt{n}}) \right)^n \\ &= \lim_{n\to\infty} \left(1 + \frac{- \frac{t^2}{2} - O(\frac{1}{\sqrt{n}})}{n} \right)^n \\ &= e^{- \frac{t^2}{2}} \end{align}$
となり、標準正規分布と特性関数が一致するので、Znは標準正規分布に収束する。

ブログ記事を読む

「統計学が最強の学問である［数学編］」のロジスティック回帰をやってみた

何年も前に買って読んでいなかった「統計学が最強の学問である」シリーズを、今更ながら読んでいる。その［数学編］はほとんど知っている内容だったので、半分くらいは流し読みしたのだが、所々、特に微分・積分以降はわかっているようでわかっていなかったことを思い知らされて勉強になったし、大変面白かった。

さて、p.490のロジスティック回帰の最急降下法の説明の所に、「気になる方はエクセルや何らかのプログラミング言語を使って挑戦していただければと思いますが、...」と書かれている。気にならない訳がないので、やってみた。

問題は、次のようなデータについて、トラブルの有無がどの程度リピート率に影響するかをロジスティック回帰分析で調べるというものである。

利用時間	トラブル	リピート	人数
ランチ	なし	なし	207
ランチ	なし	あり	23
ランチ	あり	なし	18
ランチ	あり	あり	2
ディナー	なし	なし	435
ディナー	なし	あり	290
ディナー	あり	なし	15
ディナー	あり	あり	10

●コード

import numpy as np

# 元データ(p.483)
data = np.array([
    # ディナーダミー、トラブルダミー、リピートダミー、人数
    [0, 0, 0, 207],
    [0, 0, 1, 23],
    [0, 1, 0, 18],
    [0, 1, 1, 2],
    [1, 0, 0, 435],
    [1, 0, 1, 290],
    [1, 1, 0, 15],
    [1, 1, 1, 10]])

# 各行を人数分に展開したもの
data_flat = np.repeat(data[:, :3], data[:, 3], axis=0)

X = data_flat[:, :2] # ディナーダミー、トラブルダミー
X = np.hstack((np.ones_like(X[:,:1]), X)) # 1列目に全て1の列を追加
y = data_flat[:, 2]  # リピートダミー

# コスト関数 C(β) = -log L(β) = - \sum { (y - 1) x^T β - log(1 + exp(-x^T β)) }
# L(β) = \Pi p^y (1-p)^(1-y), p = 1 / (1 + exp(-x^T β))
def cost_function(X, y, beta):
    return -np.sum((y - 1) * (X @ beta) - np.log(1 + np.exp(-X @ beta)))

# 勾配 ∂C/∂β のテスト(p.490)
beta = np.array([1.0, 1.0, 1.0])
delta = 0.0001
gradient = [
    (cost_function(X, y, beta + [delta, 0, 0]) - cost_function(X, y, beta)) / delta,
    (cost_function(X, y, beta + [0, delta, 0]) - cost_function(X, y, beta)) / delta,
    (cost_function(X, y, beta + [0, 0, delta]) - cost_function(X, y, beta)) / delta,
]

print("first gradient =", gradient)

# 勾配が最小になるβを求める
# βから勾配に学習率を掛けたものを引くのを繰り返す
learning_rate = 0.005

for i in range(120):
    gradient = np.array([
        (cost_function(X, y, beta + [delta, 0, 0]) - cost_function(X, y, beta)) / delta,
        (cost_function(X, y, beta + [0, delta, 0]) - cost_function(X, y, beta)) / delta,
        (cost_function(X, y, beta + [0, 0, delta]) - cost_function(X, y, beta)) / delta,
    ])
    beta -= learning_rate * gradient

print("beta =", beta)

# オッズ比(p.492)
print("odds of beta1, beta2 =", np.exp(beta[1:]))

# scikit-learnのLogisticRegressionで検算する
from sklearn.linear_model import LogisticRegression

X = data_flat[:, :2] # ディナーダミー、トラブルダミー
y = data_flat[:, 2]  # リピートダミー

lr = LogisticRegression().fit(X, y)
print("sklearn beta0 =", lr.intercept_, ", (beta1, beta2) =", lr.coef_[0])
print("sklearn odds of beta1, beta2 =", np.exp(lr.coef_[0]))

●実行結果

first gradient = [523.157877723861, 362.3960970844564, 29.430456188492826] beta = [-2.19686992 1.79119201 0.00347654] odds of beta1, beta2 = [5.99659623 1.00348259] sklearn beta0 = [-2.12298445] , (beta1, beta2) = [ 1.70835681 -0.00918478] sklearn odds of beta1, beta2 = [5.51988382 0.99085727]

betaの1つ目が切片、2つ目と3つ目が回帰係数で、それぞれ小数点第2位を四捨五入すると[-2.2, 1.8, 0.0]、回帰係数をオッズ比にすると[6.0, 1.0]である。p.493の結果と一致した。トラブル有無のオッズ比が1.0なので、リピート率への影響は無しという結果である。
検算の為にscikit-learnのLogisticRegressionでもやってみたら、切片が-2.1、回帰係数の1つ目が1.7と少し違う結果になり、LogisticRegressionのハイパーバラメーターを色々変えてみても本の記載と同じ結果にはできなかった。そういうものなのだろうか。コスト関数の勾配は本の結果の方が小さかったので、本の結果の方が良いと思われる。

ブログ記事を読む

吃音の治癒に関して統計解析した論文の統計用語調査

前の記事で紹介した、吃音の治癒に関して統計解析した論文
"Spontaneous" late recovery from stuttering: Dimensions of reported techniques and causal attributions
に、筆者にとって馴染みの無い統計学の用語がいくつかあったので、この機会に調べた。
この論文では主成分分析を行っているが、それらの用語は因子分析で用いられるものが多かった。

Kaiser-Meyer-Olkin (KMO) index (§2.3)

主に因子分析の文脈において、データがどれくらい因子分析に適しているかを示す指数。
データ中の各観測変数の妥当性の評価、及び分析モデルの妥当性の評価に用いられる。
MSA（Measure of Sampling Adequacy, 標本妥当性）とも言われる。

定義としては、
$KMO=\frac{\sum\sum_{i\neq j} r_{ij}^2}{\sum\sum_{i\neq j} r_{ij}^2 + \sum\sum_{i\neq j} q_{ij}^2}$
但し、
$Q=SR^{-1}S,\ S=diag(R^{-1})^{-\frac{1}{2}}$
Rはデータの観測変数の相関行列
Qは偏相関行列
r_ijはRの要素
q_ijはQの要素

KMOの3つのΣΣ部分において、iについて和を取らなければ、観測変数毎の（i番目の観測変数の）KMO indexとなる。

意味としては概ね、（相関＋偏相関）中の相関の割合であり、大きいほど観測変数間に相関がある＝共通因子があるということになる。因子分析においては0.8以上が好ましく^[3]、0.5未満は"unacceptable"とされているらしい。（0.5だと相関係数と偏相関係数が同じ大きさなので共通因子無しである）

KMO indexが低い場合、変数毎のKMO indexが低い変数を除外する方法や、データ数を増やす方法があるらしい。今回の論文では、KMO indexが0.5を超えるようにデータ数を増やすという使い方がなされた。

主に因子分析の文脈で使われるが、主成分分析でデータを縮約する場合や、共通因子を見つけようとする場合にも使われるようだ。

参考文献
[1] Kaiser-Meyer-Olkin test - Wikipedia
[2] R: Kaiser-Meyer-Olkin criterion
[3] A Modernized Heuristic Approach to Robust Exploratory Factor Analysis
[4] https://htsuda.net/stats/factor-analysis.html

polychoric correlation matrix (§2.3)

ポリコリック相関係数行列。ポリコリック相関係数は、因子分析や主成分分析で順序尺度のデータを適切に扱う為に用いられるもので、順序尺度同士の相関係数を、通常の相関係数（ピアソンの積率相関係数）とは異なり、順序尺度を等間隔ではなく、正規分布に従う連続値が何段階かに分けられたものと仮定して計算される。
今回の論文では、アンケートの結果が基本的に
"definitely not true" (1),
"probably not true" (2),
"don't know" (3),
"probably true" (4),
"definitely true" (5)
の5値の順序尺度なので、それらのポリコリック相関係数に基づいて主成分分析が行われている（ポリコリック相関係数行列の固有ベクトルを主成分としていると思われる）。

参考文献
[5] 小杉先生の資料

Parallel Analysis(PA), eigenvalues > 1 criterion (§2.3 etc.)

因子分析における因子数の決定方法。
以下は代表的な因子数の決定方法の例である。

カイザー・ガットマン基準: 相関行列の固有値の内、1以上のものの個数を因子数とする方法。
今回の論文に書かれている"eigenvalues > 1 criterion"はこれのことである。
1つの因子にしか負荷しない因子の固有値が最大1だから合理的である。
昔から統計ソフトが対応していたり、デフォルトの選択だったりしたので、よく使われていたらしい。
標本誤差があると因子数が多くなる傾向があり、あまり良い指標ではない。
スクリーテスト: 相関行列の固有値を大きい順にプロットし、下の方が成す線から離れた大きい固有値の数を因子数とする方法。
過去にはよく使われていたらしい。
平行分析(Parallel Analysis, Horn's PA): 元のデータと同じサイズの正規乱数行列の相関行列の固有値より大きい固有値の数を因子数とする方法。
現在最も推奨される方法の1つのようだ。
今回の論文では、因子数というか、データを説明するのに有効な主成分の数（次元圧縮後の主成分の数）を決めるのに用いられている。
MAP(Minimum Average Partial): 主成分の影響を取り除いた偏相関（≒誤差の相関）の二乗和が最小になる主成分の数を因子数とするような感じの方法。
最尤解のカイ二乗検定: 因子数を1から順に増やして、モデルの適合度を示すカイ二乗値が初めて有意でなくなった(p>0.05)所で止めるという方法。
因子分析のモデルを最尤法で求める場合、カイ二乗検定でモデルの適合度を調べられる。
結果がサンプルサイズによって変わるので、あまり良い指標ではないらしい。
BIC（ベイズ情報量基準）: 同じく最尤解の場合に使える方法で、BICが最小になる因子数を採用する方法。
特に順序尺度のデータをポリコリック相関係数で扱う場合に合理的らしい（ポリコリック相関係数がベイズ法で計算される為。参考文献[3]より）。

参考文献
[6] 因子分析における因子数決定法（堀先生の論文）
[7] Determining the Number of Factors to Retain in an Exploratory Factor
Analysis Using Comparison Data of Known Factorial Structure
[8] R: Kaiser-Guttman Criterion

Cohen's effect size of d (§3.1)

効果量(effect size)というのは統計的仮説検定において群間の差などを表す指標で、Cohen's dやHedges' gがよく用いられる。
Cohen's dの定義は、（群1の平均 - 群2の平均）÷プールされた標本標準偏差、であり、群間の平均の差が、帰無仮説における標準偏差の何倍かを示す。
p-valueとは異なり、サンプルサイズに依存しないので、p-valueが小さすぎてわかりにくい場合などにp-valueと共に示すのが有効である。
今回の論文でも、p-valueと共に示されている。

varimax rotation with Kaiser normalization (§3.2)

バリマックス回転(varimax rotation)は、主成分（因子）の解釈性を高める為の、多次元空間における主成分（因子）軸を回転する方法の1つ。主成分（因子）の負荷量が高い変数の個数を最小化するように、主成分（因子）負荷行列の各要素を2乗したものの各列の分散が最大になるように直交回転する。
Kaiser normalizationは、バリマックス回転において、全ての変数が回転後の解に等しく影響するように、主成分（因子）負荷を正規化する。

full information maximum likelihood estimation (§3.2)

欠損データを補完する方法の1つで、変数それぞれがある分布（正規分布など）に従っていると仮定して欠損部分を最尤推定する方法。FIMLと略されるらしい。
今回の論文では、"the correlation matrix was estimated using the full information maximum likelihood estimation."と、相関係数行列の推定にFIMLが用いられるように書かれている。例えばRのpsychパッケージには相関係数行列をFIMLで計算する関数があり^[10]、内部ではFIMLにEMアルゴリズムが用いられる^[11]ようなので、欠損データをFIMLで補完してから相関行列を求めるのではなく、欠損データの補完と相関行列の推定を同時に行ったものと推測した。

参考文献
[9] FIML Basic Concepts | Real Statistics Using Excel
[10] R: Find a Full Information Maximum Likelihood (FIML) correlation...
[11] Package 'lavaan'

Mann-Whitney U test (§3.4)

2群の中央値に差があるかどうかを検定する、ノンパラメトリックな手法。
2群のデータを全て合わせて小さい方から順位をつけ、群1,群2のそれぞれの順位の和をR₁,R₂とし、U₁,U₂を
U₁ = n₁n₂ + n₁(n₁+1)/2 - R₁
U₂ = n₁n₂ + n₂(n₂+1)/2 - R₂
（n₁,n₂はそれぞれ群1,群2のデータ数）とし、その小さい方を検定量Uとし、Uの確率分布からp-valueを求めるという感じの方法。
統計検定準1級の出題範囲にもある、お馴染みのウィルコクスンの順位和検定と同じ結果になるらしい。

Benjamini-Hochberg procedure, Bonferroni correction (§3.4)

多重検定を行うと誤検出率が上がってしまう問題に対処する方法。
同じ有意水準で複数の検定を行うと、実際には有意差が無くてももどれかでは有意となる確率("familywise error rate", FWER)が上がる。例えば有意水準α=0.05とすると、実際には有意差無しでも有意となってしまう(Type I errorの)確率が0.05なので、同じ有意水準で検定を3つ行ってどれかでは有意となる確率は 1-(1-0.05)³≒0.14 となる。

Bonferroni correction（ボンフェローニ補正）は、mを検定の数として、複数の検定における有意水準を一律にα/mとする単純な方法。m個の検定のどれかが有意となる確率はαより小さくなるので、検出の基準が厳しくなる。即ち、有意なものを検出できなくなる(Type II error)確率が上がる。mが大きいほど顕著になる。

Benjamini-Hochberg procedure（BH法）は、m個のp-valueを小さい順に並べたi番目のp-valueをp_iとして、p_i≦(i/m)αを満たす最大のiまでのp_iに対応する帰無仮説を棄却するという方法。ボンフェローニ補正の問題を解消できる。

今回の論文には、年齢の違い（中央値より大か小か）、性別の違いについてMann-WhitneyのU検定で有意差があるかどうかを調べる際に、多重検定の問題はfamilywise error rateを制御するのではなくBenjamini-Hochberg procedureでp-valueを調整することで対処している、Bonferroni correctionよりは厳しくないので偽陰性の数を下げられると書かれている。

ブログ記事を読む