最初に
取得したデータの観測値群に対して代表的な値を取るとき、「平均値と中央値はどちらがいいのか」というのが一般に良く言われると思います。
この時、「代表的な値」というのは集団を示すそれっぽい値というか、そういうのを示したいんだと思います。
例えば所得の分布を見る場合は代表値として中央値の方が良いと言われるし、以下のような釣り鐘型の分布なら平均のほうが良いよと違いを解説するページでは良くあると思うのですが、何を元にして「良い」のかを最尤法を用いて自分なりに考えてみました。
- 分布の例
推定
まず観測データを取得する母集団が何らかの確率分布に従っていると仮定します。
この時、観測データの標本は互いに独立であり、すべて同じ仮定した分布から取得されているものとします。
(左 : 母集団確率分布, 右:母集団からの観測データ)
もちろんですが、母集団確率分布はどのような形をしているか分かりません。
よって確率分布の未知のパラメータ$\theta$を推定し、その推定したパラメータ$\hat{\theta}$の値を評価したいです。
最尤法
パラメータを推定する方法の一つとして最尤法があります。
確率分布の確率密度関数を$f(x:\theta)$とした時、標本の独立同一性より、同時確率密度関数は
L(\theta) := \prod_{i=1}^{n}f(x_i:\theta)
となります。これを$\theta$の関数として扱う時尤度関数と言います。
尤度は、得られた標本が確率分布からどれだけ生成されやすいかの指標と解釈することができます。
例えば、パラメータ$\theta_1$を適当に決定し、データを10回観測してそれらがすべて400だった場合、そのパラメータ$\theta_1$の元で10回観測してそれらがすべて400になる確率はどれくらいなのか というのを見ることができます。もちろん的外れすぎるパラメータを当て嵌めた時、尤度は小さくなります。
最尤推定量
そうしたら「尤度が大きくなるパラメータを探せばいいじゃん」という風になります。
ということで、仮定した分布のパラメータ領域内において、尤度を最大化するパラメータ値をパラメータの推定量として用いることを、最尤推定量と言います。
正規分布
「平均を取るといいですよ！」代表の正規分布ですが上の例で出したような形です。
確率密度関数とパラメータは以下のようになります。
f(x | \mu, \sigma^2) = \frac{1}{\sqrt{2\pi}\sigma}exp \Bigl\{-\frac{(x - \mu)^2}{2\sigma^2}\Bigr\}
つまり母集団の確率分布を正規分布としたときの対数尤度は
L(\mu, \sigma^2) := \prod_{i=1}^{n}\frac{1}{\sqrt{2\pi}\sigma}exp \Bigl\{-\frac{(x_i - \mu)^2}{2\sigma^2}\Bigr\}
とすることができます。
総積では、計算がしづらい為、対数を取ります。これを対数尤度と言います。
\ell (\mu, \sigma^2) := \sum_{i=1}^{n} -\frac{1}{2}log(2\pi\sigma^2) -\frac{(x_i - \mu)^2}{2\sigma^2} \\
= -\frac{1}{2}n log(2\pi\sigma^2)- \sum_{i=1}^{n} \frac{(x_i - \mu)^2}{2\sigma^2}
パラメータを最大化したい場合、そのパラメータに対して偏微分し、その時の値が0となる点を求めればよいです。
$\mu$の場合
\frac{\partial}{\partial \mu}= -\frac{1}{2}n log(2\pi\sigma^2)- \sum_{i=1}^{n} \frac{(x_i - \mu)^2}{2\sigma^2} \\
\frac{\partial}{\partial \mu}=- \sum_{i=1}^{n}\frac{x_i^2 - 2x_i\mu + \mu^2}{2\sigma^2} \\
0 = - \sum_{i=1}^{n}\frac{-2x_i + 2\mu}{2\sigma^2} \\
0 = - \sum_{i=1}^{n}\frac{-x_i + \mu}{\sigma^2} \\
0 = - \sum_{i=1}^{n} -x_i + \mu \\
n\mu = \sum_{i=1}^{n} x_i \\
\mu = \frac{1}{n} \sum_{i=1}^{n} x_i
となります。
母集団確率分布が正規分布に従っているとき、観測データからパラメータ$\mu$についての尤度を最大化したい時、
\frac{1}{n} \sum_{i=1}^{n} x_i
とすればいいことが分かります。観測データから平均を取るときと一緒です。
不偏推定量
不偏性ですが、推定した値が平均的に大きすぎず小さすぎずで、母数を推定している状態の事です。
統計WEBの弓矢の図が分かりやすいです。
推定したパラメータでの期待値が母集団のパラメータと一致するよ ということです。
E[\hat{\theta}] = \theta
一方で、上記式が成り立たずずれてしまっている場合、ずれている分をバイアスと呼びます。
不偏推定量はバイアスが0の推定量という事です。
一様最小分散不偏推定量
バイアスが0の推定量が不偏推定量という事ですが推定したパラメータでの期待値が母集団のパラメータと一致すれば不偏推定量になります。
よって、分散についての評価がされていないです。
要は、以下のグラフのように、平均を取れば緑の点と一致する（バイアスが0）場合でも、ばらつきが違うときがあるよね という事です。
- 分散の違う不偏推定量(イメージ)
不偏推定量の中でもばらつきが違うものがある為、ばらつきについての評価をしたいです。要するにばらつきが一番小さい不偏推定量を探したいです。
ということで推定量のばらつき$V[\hat{\theta}]$を最小化した推定量を一様最小分散不偏推定量と言います。
クラメール・ラオの不等式
一様最小分散不偏推定量であるかどうかを判定するものとして、クラメール・ラオの不等式があります。
V[\hat{\theta}] \geq J_n(\theta)^{-1}
逆数を取っている$J_n(\theta)$はフィッシャー情報量といい
J_n(\theta) = -E\Biggl[ \frac{\partial^2}{\partial \theta^2}\ell(\theta) \Biggr]
で表すことができます。式は同時確率密度関数に対数尤度を取る部分までは最尤推定量を求める計算と同じで、推定したいパラメータについての偏微分を2回行います。
正規分布でフィッシャー情報量を求めてみます。
最尤推定量を求める際に一度$\mu$について微分したところから続けて行います。
\frac{\partial}{\partial \mu} = - \sum_{i=1}^{n}\frac{-x_i + \mu}{\sigma^2} \quad (\mu について1回微分済) \\
= - \sum_{i=1}^{n}\frac{1}{\sigma^2} \\
-E\Bigl[ - \sum_{i=1}^{n}\frac{1}{\sigma^2} \Bigr] = \frac{n}{\sigma^2}
で、正規分布のフィッシャー情報量は$\frac{n}{\sigma^2}$です。
クラメール・ラオの不等式の話に戻りますが、不偏推定量$\hat{\theta}$をどのように選んでも、その分散をフィッシャー情報量の逆数より小さくすることはできないという事です。
V[\hat{\theta}] \geq \frac{\sigma^2}{n}
$V[\hat{\theta}]$の値はいくらなの？という話になりますが、正規分布の最尤推定量は標本平均と同じでした。
よって標本平均の分散を求めればいいことになります。
標本平均の分散は
V \Bigl[ \frac{1}{n} \sum_{i=1}^{n}x_i \Bigr] = \frac{1}{n}V[x_i] = \frac{\sigma^2}{n}
となり、正規分布の場合は不偏推定量の分散$V[\hat{\theta}]$とフィッシャー情報量の逆数$J_n(\theta)^{-1}$が一致します。
このように互いの値が一致しているとき、不偏推定量$\hat{\theta}$はクラメール・ラオの下限を達成しているといい、達成している不偏推定量を有効推定量と言います。
よって、有効推定量であれば、これ以上分散の小さい推定量は存在しない為、一様最小分散不偏推定量であることが分かります。
正規分布との関係
つまりどういうこと？ってことになると思いますが、「観測データに対して平均を取る」という行為は、母集団が正規分布に従っている場合、一様最小分散不偏推定量であり、その分布の$\mu$に対して最も良い推定をしている という事になります。
よって、「釣り鐘型の分布に対しては平均を取るのが良い」のは平均が一様最小分散不偏推定量であるからです。
ちなみにですが、正規分布の場合中央値と平均値が一致するのでどっちでも良かったりします
問題なのは、裾の長い分布の場合中央値のほうが「良い」のはなんで？という事です。
裾の長い分布について
最初のページで例として出した所得の分布や、CODEの購買データで良く見る特定期間中の購買回数分布等は右裾の長い分布と呼ばれます。
- 右裾の長い分布例
多人数による特定期間の購買回数分布は負の二項分布で従うと言われており、ポアソン分布とガンマ分布を掛けた式で表されます。
導出やイメージは以下のページが分かりやすいです。
一方所得の分布は対数正規分布に従うと言われています。
裾の長い分布について期待値のパラメータを推定する際、標本平均を用いることは不偏推定量ではありますが、最小分散不偏推定量ではない時があります。
ロバストな推定量
じゃあ、「期待値のパラメータ推定に標本平均を用いた不偏推定量よりも、分散が少ない不偏推定量を手軽に探そうよ」となった時に中央値が出てきます。
ちなみにここで「最小分散不偏推定量を探そうよ」とならないのは、最尤法については仮定を置く分布によっては計算が出来なかったりするからです。
ロバストはで頑健という意味で、統計的に「一般的には，ある統計手法が仮定している条件を満たしていないときにも，ほぼ妥当な結果を与えるとき」とありました。
個人的なロバストの感覚ですけど「まぁなんか代表値にするには妥当な値」って感じです。
中央値は、平均値よりも外れ値の影響を受けづらいので、分布の裾が長く外れ値の影響を受けやすいような観測データの場合は中央値の方がより良い推定量になることがあります。
半コーシー分布
例えば、半コーシー分布のような、極端な外れ値が存在し期待値が存在しない分布からのランダムサンプリングをして分布と平均値中央値を出してみます。
| 集計指標 | 値 |
|---|---|
| 中央値 | 8.17 |
| 平均値 | 23.31 |
右側のグラフの方が分かりやすいと思いますが、中央値と平均値、どっちの方が「それっぽい値」でしょうか。
中央値の方がその値付近に観測サンプルが密集しているので、母集団分布も中央値付近の値を取りやすい と言えるのではないでしょうか。
このように、正規分布から離れた形の分布の場合、中央値の方が推定量が望ましいものになる場合があります。
ロバストな推定法いろいろ
期待値のロバストな推定法を紹介します。
中央値
言わずもがなですが中央値です。
観測データを値が小さい順に並べ替えて、真ん中の値を推定量とします。
観測データが偶数の場合は真ん中が無いので、余った2つの平均を取ります。
{Med(\,X\,) = \left\{
\begin{array}{ll}
x_{(\,\frac{n+1}{2}\,)} & ,\,nが奇数 \\
\\
\frac{x_{(\,\frac{n}{2}\,)}\,\,\,+\,\,\,x_{(\,\frac{n}{2}+1\,)}}{2} & ,\,nが偶数
\end{array}
\right.
}
最頻値
標本の中で最も観測数が多いデータの値を推定量とする方法です。
トリム平均
観測データを値が小さい順に並べ替えてから、両脇をk個除いて、その後平均を取るやり方です。
\hat{\theta}_{trim}(k) = \frac{1}{n - 2k}(x_{k + 1} + x_{k + 2} + ... + x_{n - k})
kの値はサンプル数が少なければ1か2、ある程度大きければ5%から10%程度とするのが普通だそうです。
ホッジス・レーマン推定量
観測データから2つ選んで平均を取る をすべての組み合わせに対して行います。
そうしてできた値に対して中央値を取ります。
{\mu_{HL}=Med( \{\,\frac{x_i\,+\,x_j}{2}\,\,|\,1≤i≤j≤n\,\} )
}
正規分布にしたがっているのかを判断する
最後に「正規分布に従っているのかどうか」を説明します。
結局のところ、「代表的な値」を求めるときに、正規分布なら平均、裾が長くて正規分布っぽくないなら中央値とかロバストな推定量でというのは分かったので、「正規分布かどうか」が知りたいところです。
シャピロウィルク検定やQQプロットでの可視化確認等も有りますが、歪度と尖度を用いた判断方法について書きます。
観測値$X_1,...,X_n$に基づいて、仮説を
H0 : 観測値$X_1,...,X_n$は正規分布に従う
H1 : 観測値$X_1,...,X_n$は正規分布に従わない
とし、観測データの歪度と尖度を求めます。
歪度と尖度
歪度は正規分布を0とし、片側の裾が長いと値が変動していきます。
尖度は正規分布を3とし、両側の裾の短長で値が変動していきます。（尖度は漢字は「尖る」ですが、見た目が尖れば尖るほど値が高くなるわけではないです）
図としてはこれも統計WEBが分かりやすいです。
導出方法としては、平均周りの$r$次モーメントを求めればいいのですが、詳しくは長くなるので割愛します。
正規分布の場合の積率母関数（モーメント母関数）は以下の式で表されます。
M(t) = E[e^{tX}] \\
= \int_{-\infty}^{\infty}e^{tx}f(x) dx \\
= \int_{-\infty}^{\infty}\frac{1}{\sqrt{2\pi}\sigma}exp \Bigl\{-\frac{(x - \mu)^2}{2\sigma^2} + tx\Bigr\}\\
上記式ですが、$M(t)'$とし、$t=0$を代入すると$E[X]$となります。
$M(t)''$とし、$t=0$を代入すると$E[X^2]$となります。
$M(t)^{(r)}$とし、$t=0$を代入すると$E[X^r]$となります。
歪度は3次モーメントの値を標準偏差の3乗で割った値になります。
尖度は4次モーメントの値を標準偏差の4乗で割った値になります。
- 歪度
\frac{E[X^3]}{V[X]^{\frac{3}{2}}}
- 尖度
\frac{E[X^4]}{V[X]^{2}}
分母の$V[X]$ですが、1次モーメントと2次モーメントを利用すれば求まります。
分散を求める式は
V[X] = E[X^2] - \Bigl(E[X] \Bigr)^2
となり、2次モーメント - 1次モーメントの2乗で求まります。
正規分布$V[X]$は$\sigma^2$である為、分母は以上のようになります。
歪度と尖度と求め、それらの絶対値が0(尖度は-3し正規分布の尖度が0になる様補正する)よりも大きく離れているようであれば、帰無仮説H0を棄却すれば良いです。
まとめ
- 正規分布っぽければ平均が「良い」
- 右裾・左裾が長い分布は中央値や、他のロバストな推定量を用いると平均よりも「良い」代表値がとれる
かなと思います。