Python初心者の備忘録 #13 ~統計学入門編03~

# はじめに
今回私は最近はやりのchatGPTに興味を持ち、深層学習について学んでみたいと思い立ちました！  
深層学習といえばPythonということなので、最終的にはPythonを使って深層学習ができるとこまでコツコツと学習していくことにしました。  
ただ、勉強するだけではなく少しでもアウトプットをしようということで、備忘録として学習した内容をまとめていこうと思います。  
この記事が少しでも誰かの糧になることを願っております！
**※投稿主の環境はWindowsなのでMacの方は多少違う部分が出てくると思いますが、ご了承ください。**
最初の記事：[Python初心者の備忘録 #01](https://qiita.com/Yushin-Tati/items/961dc47d6163f944f7e9)
前の記事：[Python初心者の備忘録 #12 ~統計学入門編02~](https://qiita.com/Yushin-Tati/items/af085499c49d244068ee)
次の記事：[Python初心者の備忘録 #14 ~統計学入門編04~](https://qiita.com/Yushin-Tati/items/8aeb2c190e3494dd6b78)

本記事は**推測統計の前提知識、区間推定**についてまとめてあります。
## ■学習に使用している資料
Udemy：[米国データサイエンティストが教える統計学超入門講座【Pythonで実践】
](https://www.udemy.com/course/python-stats/)
## ■推測統計入門
### ▶推測統計とは
- 推測統計には大きく「推定」と「検定」がある
推定：標本の特徴から母集団の特徴を推測する
検定：「標本同士の特徴を比較し、同じ母集団から抽出したものなのか」ということや「標本の特徴から母集団は正規分布かどうか」といったことをテストする
- 推定には「点推定」と「区間推定」がある
点推定：バシッと1つの値で母集団の特徴を推定(推定ではほとんど使用しない)
区間推定：ある程度の幅を持って、母集団の特徴を推定(母数は○○～△△の間かな？)
### ▶母集団(population)と標本(sample)
- 母集団(日本の成人男性)の特徴(平均年収)を標本(ランダムに抽出した成人男性1000人)から推測するような作業を「推測統計」という
- 完全な無作為抽出は困難なので、できるだけランダムになるように努力する
- すでに標本が存在する場合は、無理のない範囲の母集団を設定する必要がある
(例：1日のチップの平均 -> 「×：そのレストランの1日のチップの平均」「〇：ある曜日のチップの平均」)
- 標本から母集団の推測が主な作業となるので、母集団の設定は慎重に行う必要がある
### ▶標本統計量と母数
- 分布の特徴を表す値を統計量という(平均、分散など)
- 標本から計算された統計量を標本統計量、母集団の統計量を母数という

![stats_population_sample.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3448232/04216cc7-8728-8d20-725d-48cbc2b2fb79.jpeg)
### ▶標本分布(sampling distribution)
- 標本の中身は確率的に変動する(抽出するたびに違う標本が得られる)
- 標本統計量も同様の理由で確率的に変動する(つまり確率変数)
- 標本統計量が従う確率分布を標本分布という

例：全国の成人男性を無作為抽出したときの各標本の身長の平均
![stats_sampling_distribution.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3448232/bae07bb1-db25-9269-991e-60d032121927.jpeg)
### ▶推定量(estimator)
- 母数の推定に使用する標本統計量のこと
- 実際の値を推定値という
(標本平均：母数平均の推定量、173㎝：推定値)
- 標本分散と標本標準偏差は母数の推定量に使用できない

![stats_estimator.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3448232/9bc44ef8-76e2-27c3-41b6-f4c325f4f540.jpeg)
### ▶不偏性(unbiasedness)
- 推定量が**平均的に**母数と一致する場合、その推定量は「不偏性がある」といい、不偏推定量(unbiased estimator)と呼ぶ
- 確率を扱う文脈では平均のことを期待値という
※不偏推定量は平均、比率、不偏分散が該当する

![stats_unbiasedness.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3448232/b981b461-26ad-1b10-0003-1d82645c287a.jpeg)
母集団の平均が「171cm」だった時、標本平均も大量にとれば推定値の平均が「171cm」になることが感覚的にわかる。
つまり、**推定量が「平均的」に母数と一致する**といえるので、標本平均は不偏推定量である。
### ▶平均の標本分布
- 平均の標本分布の平均は母平均と同じμになり、分散は母分散$σ^2$をnで割った値($\frac{σ^2}{n}$)、標準偏差は$\sqrt{\frac{σ^2}{n}}$となる

![stats_sample_var_std.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3448232/d08451d1-a922-1e71-b210-60db08ec8b24.jpeg)
1回に取る標本数nが大きくなればなるほど、母集団の総数に近づくので、標本平均が母平均に一致する確率が大きくなる。
つまり、分散や標準偏差の値が小さくなるということ。
#### Challenge
データセットを母集団と見立てて、50個の標本を100回抽出したときの標本平均分布(平均と分散)を確認する。
- 今回は`'tip'`について
- `df.sample(n)`でdfからランダムにn行抽出することができる
<details><summary>解答例</summary>

```python
m = 100
n = 50
column = 'tip'
sample_means = []
for i in range(m):
    sampled_df = df.sample(n)
    sample_mean = sampled_df[column].mean()
    sample_means.append(sample_mean)

#平均
np.mean(sample_means)

# 母平均
df[column].mean()

# 分散 (実際には不偏分散stats.tvar()を使った方が良い)
# 外れ値の影響で母分散と差が出てしまう
np.var(sample_means)

# 母分散/n
np.var(df[column])/n
```
</details>

### ▶不偏分散(unbiased variance)
- 母分散の不偏推定量
- 標本不偏分散の標本分布の平均は母分散$σ^2$と一致する
- 分散の式でnではなくn-1で割った値で、$s'^2$と表す

分散の式：$s^2=\frac{1}{n}\sum_{i=1}^n(x_i-\bar{x})^2$

不偏分散の式：$s'^2=\frac{1}{n-1}\sum_{i=1}^n(x_i-\bar{x})^2$
#### なぜn-1で割るのか？
1. 分散の計算の際に$\sum_{i=1}^n(x_i-\bar{x})^2$を求めているが、**この値(平方和)は$\bar{x}$が平均値の時最小となる**。
1. 本来であれば母平均μを使って分散を出したほうがいいのだが、わからないので標本平均$\bar{x}$を使用している。
1. 「1.」の特性上、分散を計算する際にμを使用したときよりも$\bar{x}$を使用したときのほうが値が小さくなる。(**標本分散が母分散より必ず小さくなる** -> 過小評価される)
$\frac{1}{n}\sum_{i=1}^n(x_i-\bar{x})^2\leq\frac{1}{n}\sum_{i=1}^n(x_i-μ)^2$
1. そのため、nではなく少し小さい値の**n-1で割ることで、過小評価の補正を行っている**。

![stats_unbiased_variance.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3448232/06b0bfb9-d2fa-3b79-3f5a-b9bc26f73156.jpeg)
### ▶不偏分散の平方根について
- 不偏分散の平方根は母集団の標準偏差σの不偏推定量にはならない
- 標準偏差(不偏分散の平方根)s'の標本分布の平均は母標準偏差σにはならない
※不偏推定量にはならないが、**推定量としては十分使用できる**。
```python
#不偏分散(下記のどちらでもOK)
stats.tvar(df['tip'])  # -> 1.914454638062471
np.var(df['tip'], ddof=1)  # -> 1.9144546380624725

# 不偏分散の平方根(下記のどれでもでもOK)
np.sqrt(stats.tvar(df['tip']))  # -> 1.3836381890011822
stats.tstd(df['tip'])  # -> 1.3836381890011822
np.std(df['tip'], ddof=1)  # -> 1.3836381890011826
```
#### Challenge
不偏分散の不偏性を実際に確認する
- データセットを母集団と見立てて、50個の標本を100回抽出して標本の分散と不偏分散の平均がそれぞれどうなるか
- 今回は`'tip'`を使用
<details><summary>解答例</summary>

```python
n = 50
m = 1000
column = 'tip'
sample_vars = []
sample_uvars = []
for i in range(m):
    sample_df = df[cplumn].sample(n)
    sample_vars.append(sample_var)
    sample_uvars.append(sample_uvar)

# 母分散
np.var(df[column])  # -> 1.9066085124966428

# 標本の不偏分散の平均
np.mean(sample_uvars)  # -> 1.9031641021632655

# 標本の分散の平均
np.mean(sample_vars)  # -> 1.8651008201200001

# 不偏分散の平均のほうが母分散に近い
```
</details>

## ■区間推定
### ▶比率の区間推定
- 標本の比率の値から母集団の比率の値を推定する
- 実際に比率の区間推定が応用されている例は多い(逝去速報、内閣支持率、視聴率など)

例：世論の内閣支持率
![stats_percent_estimation.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3448232/9d4dd3e4-d9f7-03ca-0193-1d6c937309d8.jpeg)
- 母数の値を確度をつけた区間で推定する(95%の確度で母比率は0.5~0.7)
- 上記の区間を信頼区間(CI：confidence interval)という(95%信頼区間で母比率は0.5~0.7)
※何度も標本を抽出した際に100回中5回の確率で区間の指定を間違えるということ
### ▶比率の標本分布
- 標本推定量の標本分布を使用する(例：比率の推定であれば標本比率)
※標本比率がどれほどの確率で得られるかを考える
- 比率の標本分布は二項分布から考えることができる

![stats_interval_binomial.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3448232/72e06ef1-d09f-fcfb-b7a5-f461f3050aca.jpeg)
- 二項分布はnが十分に大きいとき(2,30でOK)、正規分布に近似可能(この時平均はnp、分散はnpq)
- 比率の標本分布はnが十分大きいとき、平均p、分散$\frac{pq}{n}$の正規分布に近似できる
※母比率pと標本分布の平均pが一致するので、標本比率は母比率の不偏推定量といえる

![stats_sample_interval_binomial.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3448232/47500833-1254-1aa1-74be-8269dfc92062.jpeg)
### ▶比率の標本分布の標準化
- 現状では区間推定の信頼区間などの計算が面倒なので、標準化して扱いやすくする
- 標準化の式は$z=\frac{x-\bar{x}}{s}$なので、比率の場合$z=\frac{x/n-p}{\sqrt{pq/n}}$となる
- この時、95%信頼区間を求める場合は、zが-1.96～1.96の区間と考えればいい

![stats_interval_standard.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3448232/9909cdd3-aa4b-18e2-bfbd-be98407246db.jpeg)
`stats.binom.interval(alpha, n, p)`で比率の区間推定を行うことができる。
`alpha`は信頼区間、`n`は標本数、`p`は標本の比率を指定する。
※返り値は比率ではなく、数の区間なので注意
```python
# 比率の区間推定
stats.binom.interval(0.95, 1000, 0.6)  # -> (570.0, 630.0)

# つまり、母比率は95%信頼区間で57%~63%の間にあるということになる
```
#### Challenge
実データで100回比率の区間推定を行い、母比率を何回正しく推定できているか確認する。
- 95%信頼区間
- `'sex'`を使用
<details><summary>解答例</summary>

```python
n = 50
m = 100
alpha = 0.95
# 母数の男性比
population_male_ratio = df['sex'].value_counts(normalize=True)['Male']
results = []
for _ in range(m):
    sample_df = df.sample(n)
    male_ratio = sample_df.value_counts(normalize=True)['Male']
    # 区間推定
    min_num, max_num = stats.binom.interval(alpha, n=n, p=male_ratio)
    min_ratio = min_num / n
    max_ratio = max_num / n
    # True, Falseの条件作成
    correct = min_ratio < population_male_ratio < max_ratio
    results.append({'min': min_ratio, 'max': max_ratio, 'correct':correct})

# ListをDataFlameへ
result_df = pd.DataFrame(results)
# True, Flaseの数を取得
result_df['correct'].value_counts()
```
</details>

### ▶平均値の区間推定
- 基本的な考え方は比率の区間推定と同じで、平均の標本分布について考える
- 母集団が正規分布だった場合は、平均の標本分布も正規分布となる
※[平均の標本分布](#▶平均の標本分布)参照
- 母集団が正規分布じゃない場合でも、標本数nが十分に大きい(5個とかでもいいのであまり気にしなくていい)場合は平均の標本分布は正規分布となる
※**中心極限定理**によるもの
### ▶平均標本分布の標準化
- 比率同様、扱いやすくするために標準化を行う
- $z=\frac{x-\bar{x}}{s}$より、$z=\frac{\bar{x}-μ}{\frac{σ}{\sqrt{n}}}$で標準化を行う

![stats_mean_interval_normalize.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3448232/ff1134bf-e042-b774-69d1-0104fa6c4d46.jpeg)
- 母標準偏差σはわからないので、不偏分散の平方根s'で代用する
※s'はあくまで推定量なので標準正規分布は本来のものからずれた分布(t分布と呼ばれる)になるが、nが十分大きければ近似できる(n=30くらい)
- zが-1.96～1.96の間で95%信頼区間となる

![stats_mean_interval_normalize_formula.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3448232/b4e1e1c2-8707-49c8-3c60-5ccef5e0c88a.jpeg)
`stats.norm.interval(alpha, loc, scale)`で平均の区間推定が可能。
`alpha`は信頼区間、`loc`は標本平均、`scale`は$\frac{s'}{\sqrt{n}}$を指定する。
```python
# 平均の区間推定
sample_df = df['total_bill'].sample(100)
ustd = stats.tstd(sample_df)
loc = np.mean(sample_df)
scale = ustd / np.sqrt(100)
stats.norm.interval(0.95, loc=loc, scale=scale)  # -> (18.036556186791557, 21.668843813208426)
```
#### Challenge
50個の標本平均から母平均の区間推定を95%信頼区間で100回行い、ない位階正しく推定できるか確認する。
- 今回は`'tip'`を使用
<details><summary>解答例</summary>

```python
column = 'tip'
alpha = 0.95
n = 50
m = 100
population_mean = df[column].mean()
results = []
for _ in range(m):
    sample_df = df.sample(n)
    sample_mean = sample_df[column].mean()
    sample_var = stats.tvar(sample_df[column])
    min_val, max_val = stats.norm.interval(alpha, loc=sample_mean, scale=np.sqrt(sample_var/n))
    correct = min_val < population_mean < max_val
    results.append({'min': min_val, 'max': max_val, 'correct': correct})

result_df = pd.DataFrame(results)
result_df['correct'].value_counts()
```
</details>

### ▶t分布
- $z=\frac{\bar{x}-μ}{\frac{σ}{\sqrt{n}}}$で標準化した標本分布は標準正規分布(z分布)となるがσの代わりに不偏分散の平方根s'を使用して標準化した標本分布はt分布と呼ばれる
- 95%信頼区間も-1.96～1.96ではなくなる
- 自由度(df：degree of freedom)n-1をパラメータに取り、標準正規分布よりも裾野が広くなる

![stats_t_df.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3448232/ed0307a3-a284-8f34-a16d-6a00dc984cd9.jpeg)
- nが十分に大きければ、標準正規分布に近似することができる(大体n≧30ぐらい)

![stats_t_df_sample.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3448232/1647d74e-887a-5c9b-a73e-d6a77c996e90.jpeg)
`stats.t(df)`でt分布の値を出力することができる。(`df`はDataFrameではなく、自由度のことなので注意)
`stats.t(df).pdf(x)`などで値xの際のt分布の値を出力することができる。
```python
# 自由度1、6、11のt分布と標準正規分布を描画
x = np.linspace(-3, 3, 100)
z = stats.norm.pdf(x)
for dof in range(1, 12, 5):
    t = stats.t(dof).pdf(x)
    plt.plot(x, t, label=f"t dist(dof={dof})")
plt.plot(x, z, label="std norm", linewidth=3)
plt.legend()
```
![stats_t_and_std.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3448232/1f079d68-4211-9198-7bb2-d01304c7229b.jpeg)
#### Challenge
標本数50の標本平均から正規分布と使った母平均の区間推定とt分布を使った母平均の区間推定で結果がどう変化するか確認する
- 95%信頼区間
- `'tip'`を使用
- stats.t.interval(alpha, loc, scale, df)でt分布の区間推定が可能
<details><summary>解答例</summary>

```python
column = 'tip'
alpha = 0.95
n = 50
population_mean = df[column].mean()

sample_df = df.sample(n)
sample_mean = sample_df[column].mean()
sample_var = stats.tvar(sample_df[column])
min_val_norm, max_val_norm = stats.norm.interval(alpha, loc=sample_mean, scale=np.sqrt(sample_var/n))
min_val_t, max_val_t = stats.t.interval(alpha, loc=sample_mean, scale=np.sqrt(sample_var/n), df=n-1)

# 正規分布の結果
print(min_val_norm, max_val_norm)

# t分布の結果
print(min_val_t, max_val_t)
```
</details>

#### 次の記事
[Python初心者の備忘録 #14 ~統計学入門編04~](https://qiita.com/Yushin-Tati/items/8aeb2c190e3494dd6b78)