平均値と中央値の使いどころ

# 最初に
取得したデータの観測値群に対して代表的な値を取るとき、「平均値と中央値はどちらがいいのか」というのが一般に良く言われると思います。
この時、「代表的な値」というのは集団を示すそれっぽい値というか、そういうのを示したいんだと思います。

例えば[所得の分布](https://www.mhlw.go.jp/toukei/saikin/hw/k-tyosa/k-tyosa09/2-2.html)を見る場合は代表値として中央値の方が良いと言われるし、以下のような釣り鐘型の分布なら平均のほうが良いよと違いを解説するページでは良くあると思うのですが、何を元にして「良い」のかを最尤法を用いて自分なりに考えてみました。

- 分布の例

![分布例.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/635079/68d1c917-38b0-cfcc-ec0d-7599f95f084f.png)


# 推定
まず観測データを取得する母集団が何らかの確率分布に従っていると仮定します。
この時、観測データの標本は互いに独立であり、すべて同じ仮定した分布から取得されているものとします。

(左 : 母集団確率分布, 右:母集団からの観測データ)
![母集団分布例.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/635079/bccf68e8-7ee9-4269-571a-4ae4852feee1.png)


もちろんですが、母集団確率分布はどのような形をしているか分かりません。
よって確率分布の未知のパラメータ$\theta$を推定し、その推定したパラメータ$\hat{\theta}$の値を**評価**したいです。



## 最尤法

パラメータを推定する方法の一つとして最尤法があります。
確率分布の確率密度関数を$f(x:\theta)$とした時、標本の独立同一性より、同時確率密度関数は


```math
L(\theta) := \prod_{i=1}^{n}f(x_i:\theta)
```
となります。これを$\theta$の関数として扱う時尤度関数と言います。
尤度は、得られた標本が確率分布からどれだけ生成されやすいかの指標と解釈することができます。
例えば、パラメータ$\theta_1$を適当に決定し、データを10回観測してそれらがすべて400だった場合、そのパラメータ$\theta_1$の元で10回観測してそれらがすべて400になる確率はどれくらいなのか というのを見ることができます。もちろん的外れすぎるパラメータを当て嵌めた時、尤度は小さくなります。


## 最尤推定量

そうしたら「尤度が大きくなるパラメータを探せばいいじゃん」という風になります。
ということで、仮定した分布のパラメータ領域内において、尤度を最大化するパラメータ値をパラメータの推定量として用いることを、**最尤推定量**と言います。

### 正規分布

「平均を取るといいですよ！」代表の正規分布ですが上の例で出したような形です。
<img width="470" alt="標準正規分布.png (23.5 kB)" src="https://img.esa.io/uploads/production/attachments/10574/2021/09/16/62466/d39a918b-4a8a-45a7-8669-403735c07a8e.png">

確率密度関数とパラメータは以下のようになります。

```math
f(x | \mu, \sigma^2) = \frac{1}{\sqrt{2\pi}\sigma}exp \Bigl\{-\frac{(x - \mu)^2}{2\sigma^2}\Bigr\}
```

つまり母集団の確率分布を正規分布としたときの対数尤度は

```math
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\}
```

とすることができます。
総積では、計算がしづらい為、対数を取ります。これを対数尤度と言います。

```math
\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$の場合

```math
\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$についての尤度を最大化したい時、

```math
\frac{1}{n} \sum_{i=1}^{n} x_i
```
とすればいいことが分かります。観測データから**平均を取るときと一緒**です。

## 不偏推定量

**不偏性**ですが、推定した値が平均的に大きすぎず小さすぎずで、母数を推定している状態の事です。
[統計WEBの弓矢の図](https://bellcurve.jp/statistics/course/8612.html)が分かりやすいです。
推定したパラメータでの期待値が母集団のパラメータと一致するよ ということです。


```math
E[\hat{\theta}] = \theta
```

一方で、上記式が成り立たずずれてしまっている場合、ずれている分を**バイアス**と呼びます。
不偏推定量はバイアスが0の推定量という事です。

### 一様最小分散不偏推定量

バイアスが0の推定量が**不偏推定量**という事ですが推定したパラメータでの**期待値**が母集団のパラメータと一致すれば不偏推定量になります。
よって、分散についての評価がされていないです。
要は、以下のグラフのように、平均を取れば緑の点と一致する（バイアスが0）場合でも、ばらつきが違うときがあるよね という事です。

- 分散の違う不偏推定量(イメージ)
![どちらも不偏推定量の例.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/635079/f9632006-5b1f-d74f-d75f-b1475fd6aa2d.png)


不偏推定量の中でもばらつきが違うものがある為、ばらつきについての評価をしたいです。要するにばらつきが一番小さい不偏推定量を探したいです。
ということで推定量のばらつき$V[\hat{\theta}]$を最小化した推定量を**一様最小分散不偏推定量**と言います。

## クラメール・ラオの不等式

一様最小分散不偏推定量であるかどうかを判定するものとして、クラメール・ラオの不等式があります。

```math
V[\hat{\theta}] \geq J_n(\theta)^{-1}
```

逆数を取っている$J_n(\theta)$は**フィッシャー情報量**といい

```math
J_n(\theta) = -E\Biggl[ \frac{\partial^2}{\partial \theta^2}\ell(\theta) \Biggr]
```

で表すことができます。式は同時確率密度関数に対数尤度を取る部分までは最尤推定量を求める計算と同じで、推定したいパラメータについての偏微分を2回行います。
正規分布でフィッシャー情報量を求めてみます。
最尤推定量を求める際に一度[$\mu$について微分したところ](#$\mu$の場合)から続けて行います。

```math
\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}$をどのように選んでも、**その分散をフィッシャー情報量の逆数より小さくすることはできない**という事です。

```math
V[\hat{\theta}] \geq \frac{\sigma^2}{n}
```

$V[\hat{\theta}]$の値はいくらなの？という話になりますが、正規分布の最尤推定量は[標本平均](#$\mu$の場合)と同じでした。
よって標本平均の分散を求めればいいことになります。
標本平均の分散は

```math
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の購買データで良く見る特定期間中の購買回数分布等は右裾の長い分布と呼ばれます。

- 右裾の長い分布例
<img width="600" alt="分布例_2.png (29.1 kB)" src="https://img.esa.io/uploads/production/attachments/10574/2021/09/17/62466/b779e044-f74b-4dee-b884-f1d77748b532.png">



多人数による特定期間の購買回数分布は負の二項分布で従うと言われており、ポアソン分布とガンマ分布を掛けた式で表されます。
導出やイメージは以下のページが分かりやすいです。

- [購入頻度が負の二項分布に従うと考えられているのはなぜか](https://elsur.jpn.org/diary/2020/05/24/814/)
- [可視化で理解する「負の二項分布」](https://hoxo-m.hatenablog.com/entry/20151012/p1)

一方所得の分布は対数正規分布に従うと言われています。

- [対数正規分布の仕組み](https://rikunora.hatenablog.com/entry/20100418/p1)

裾の長い分布について期待値のパラメータを推定する際、標本平均を用いることは不偏推定量ではありますが、**最小分散不偏推定量ではない**時があります。

## ロバストな推定量
じゃあ、「期待値のパラメータ推定に標本平均を用いた不偏推定量よりも、分散が少ない不偏推定量を手軽に探そうよ」となった時に**中央値**が出てきます。
ちなみにここで「最小分散不偏推定量を探そうよ」とならないのは、最尤法については仮定を置く分布によっては計算が出来なかったりするからです。

**ロバスト**はで頑健という意味で、統計的に「一**般的には，ある統計手法が仮定している条件を満たしていないときにも，ほぼ妥当な結果を与えるとき**」と[ありました](https://www.weblio.jp/content/%E9%A0%91%E5%81%A5%E6%80%A7)。
個人的なロバストの感覚ですけど「まぁなんか代表値にするには妥当な値」って感じです。

中央値は、平均値よりも外れ値の影響を受けづらいので、分布の裾が長く外れ値の影響を受けやすいような観測データの場合は中央値の方がより良い推定量になることがあります。

### 半コーシー分布

例えば、[半コーシー分布](https://ja.wikipedia.org/wiki/%E3%82%B3%E3%83%BC%E3%82%B7%E3%83%BC%E5%88%86%E5%B8%83)のような、極端な外れ値が存在し期待値が存在しない分布からのランダムサンプリングをして分布と平均値中央値を出してみます。

<img width="846" alt="半コーシー分布ランダムサンプル.png (23.1 kB)" src="https://img.esa.io/uploads/production/attachments/10574/2021/09/17/62466/76b49444-4bba-49fb-98e6-44fca2d14787.png">

| 集計指標| 値 |
| --- | --- |
| 中央値 | 8.17 |
| 平均値 | 23.31 |

右側のグラフの方が分かりやすいと思いますが、中央値と平均値、どっちの方が「それっぽい値」でしょうか。
中央値の方がその値付近に観測サンプルが密集しているので、母集団分布も中央値付近の値を取りやすい と言えるのではないでしょうか。
このように、**正規分布から離れた形の分布の場合、中央値の方が推定量が望ましいものになる場合があり**ます。

## ロバストな推定法いろいろ
期待値のロバストな推定法を紹介します。

### 中央値
言わずもがなですが中央値です。
観測データを値が小さい順に並べ替えて、真ん中の値を推定量とします。
観測データが偶数の場合は真ん中が無いので、余った2つの平均を取ります。

```math
{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個除いて、その後平均を取るやり方です。

```math
\hat{\theta}_{trim}(k) = \frac{1}{n - 2k}(x_{k + 1} + x_{k + 2} + ... + x_{n - k})
```

kの値はサンプル数が少なければ1か2、ある程度大きければ5%から10%程度とするのが普通だそうです。



### ホッジス・レーマン推定量
観測データから2つ選んで平均を取る をすべての組み合わせに対して行います。
そうしてできた値に対して中央値を取ります。

```math
{\mu_{HL}=Med( \{\,\frac{x_i\,+\,x_j}{2}\,\,|\,1≤i≤j≤n\,\} )



}
```



# 正規分布にしたがっているのかを判断する
最後に「正規分布に従っているのかどうか」を説明します。
結局のところ、「代表的な値」を求めるときに、正規分布なら平均、裾が長くて正規分布っぽくないなら中央値とかロバストな推定量でというのは分かったので、「正規分布かどうか」が知りたいところです。
[シャピロウィルク検定](https://ja.wikipedia.org/wiki/%E3%82%B7%E3%83%A3%E3%83%94%E3%83%AD%E2%80%93%E3%82%A6%E3%82%A3%E3%83%AB%E3%82%AF%E6%A4%9C%E5%AE%9A)や[QQプロット](https://bellcurve.jp/statistics/glossary/2071.html)での可視化確認等も有りますが、歪度と尖度を用いた判断方法について書きます。

観測値$X_1,...,X_n$に基づいて、仮説を

H0 : 観測値$X_1,...,X_n$は正規分布に従う
H1 : 観測値$X_1,...,X_n$は正規分布に従わない

とし、観測データの歪度と尖度を求めます。

## 歪度と尖度

歪度は正規分布を0とし、片側の裾が長いと値が変動していきます。
尖度は正規分布を3とし、両側の裾の短長で値が変動していきます。（尖度は漢字は「尖る」ですが、見た目が尖れば尖るほど値が高くなるわけではないです）
図としてはこれも[統計WEBが分かりやすい](https://bellcurve.jp/statistics/course/17950.html)です。


導出方法としては、平均周りの$r$次モーメントを求めればいいのですが、詳しくは長くなるので割愛します。
正規分布の場合の積率母関数（モーメント母関数）は以下の式で表されます。


```math
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乗で割った値になります。

-  歪度

```math
\frac{E[X^3]}{V[X]^{\frac{3}{2}}}
```

- 尖度

```math
\frac{E[X^4]}{V[X]^{2}}
```

分母の$V[X]$ですが、1次モーメントと2次モーメントを利用すれば求まります。
分散を求める式は

```math
V[X] = E[X^2] - \Bigl(E[X] \Bigr)^2
```

となり、2次モーメント - 1次モーメントの2乗で求まります。
正規分布$V[X]$は$\sigma^2$である為、分母は以上のようになります。
歪度と尖度と求め、それらの絶対値が0(尖度は-3し正規分布の尖度が0になる様補正する)よりも大きく離れているようであれば、帰無仮説H0を棄却すれば良いです。




# まとめ

- 正規分布っぽければ平均が「良い」
- 右裾・左裾が長い分布は中央値や、他のロバストな推定量を用いると平均よりも「良い」代表値がとれる

かなと思います。


## 参考資料
- [自然科学の統計学](http://www.utp.or.jp/book/b300859.html)
- [現代数理統計学の基礎](https://www.kyoritsu-pub.co.jp/bookdetail/9784320111660)
- [統計学実践ワークブック](https://www.gakujutsu.co.jp/product/978-4-7806-0852-6/)