Pythonでちゃんと統計・確率の基礎をやってみた（確率編1/3）

# はじめに
以下の書籍を参考にしながら，統計の基礎を学び直した軌跡，確率編1/3です．

https://www.amazon.co.jp/Python%E3%81%A7%E7%90%86%E8%A7%A3%E3%81%99%E3%82%8B%E7%B5%B1%E8%A8%88%E8%A7%A3%E6%9E%90%E3%81%AE%E5%9F%BA%E7%A4%8E-%E8%B0%B7%E5%90%88-%E5%BB%A3%E7%B4%80-ebook/dp/B07HHV6LL8/ref=tmm_kin_swatch_0?_encoding=UTF8&qid=1657859730&sr=8-1

この記事は5部構成の2つ目になっています．以下も参照ください．
- [準備〜データ編](https://qiita.com/o93/items/91feef63ea70070462d9)・・・投稿済み
- **確率編1**・・・この記事
- [確率編2](https://qiita.com/o93/items/a2e02cf98126f47cd66c)・・・投稿済み
- 確率編3・・・まだ
- 推測編・・・まだ

この記事を書いた動機，注意点，使用している技術，記事の構成，準備については[準備〜データ編](https://qiita.com/drafts/91feef63ea70070462d9)に掲載されています．

これは主観でしかないのですが，エンジニアにも馴染みのある「データ」と比べて「確率」は扱われる機会が少なく，理解のハードルが高いかと思います．引き続き書籍を参考にしながら，理解を助けるために補足したり，独断で省略している部分もありますので，書籍に沿って進めたい方は書籍を購入ください！

また，最初は「確率編」として3部構成を考えていたのですが，確率の基本部分を厚くしたため，以下のように確率編を3つに分けることにしました．

- **確率編1/3**
    - 推測統計
    - 確率モデル
    - 離散型確率変数
        - 1次元
        - 2次元
- 確率編２/3
    - 代表的な離散型確率分布
    - 連続型確率変数
        - 1次元
        - 2次元
- 確率編3/3
    - 代表的な連続型確率分布
    - 多次元確率変数

# 基本
前記事の[構成](https://qiita.com/o93/items/91feef63ea70070462d9#%E6%A7%8B%E6%88%90)にて，統計解析のスキルを得るには確率の知識が必要という点に触れましたが，まず，統計解析の全体像から各用語を交えながら掘り下げることで，何故確率の知識が必要なのかを具体的に説明していきます．

## 記述統計と推測統計
手元にあるデータの統計的性質を記述する枠組みを**記述統計**と言います．[データ編](https://qiita.com/drafts/91feef63ea70070462d9)に掲載した内容が，まさに記述統計です．しかし，すべてのデータを手元に揃えて分析できることは稀です．

例えば，Webサイトのアクセス履歴を過去1ヶ月分手元に得たとしても，それはWebサイトを立ち上げてから1ヶ月前までのアクセス履歴や，現在からWebサイトが閉鎖されるまでの未来のデータを含んでいません．仮に我々が過去3ヶ月〜未来3ヶ月の半年間のWebサイトの統計的性質を知りたい場合，手元にある過去1ヶ月分のデータだけを使い，まだ確定していない未来を含めた半年分の統計的性質を推測することになります．

このように，一部のデータから必要なデータ全体の統計的性質を推測する枠組みを**推測統計**と言い，今回の記事で扱うことにします．

## サンプルデータ
上記の例のようなデータを扱うのは複雑なので，書籍と同様に，以下のようなデータを用いて推測統計を説明していきます．
- 全体のデータ・・・高校の全生徒`400`人の数学のテスト結果
- 手元のデータ・・・生徒Aさんが聞いて回った`20`人分のテスト結果

Aさんは高校の生徒で，全生徒のテスト結果を入手できないという制限があります．そんな条件下で，自分の点数の良し悪しを知るため全体の平均点を推測したい！というモチベーション（動機・目的）があります．図にすると以下のような感じ．

![スクリーンショット 2022-07-24 16.42.48.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/44273/5576e6bd-2438-8bb1-6dbc-483ad3080319.png)

もちろん実際に生徒に聞いて回るわけにはいかないので，サンプルデータとプログラムを用います，，，

## 推測統計の全体像
先程挙げた図が推測統計の全体像になります．が，推測統計を説明する上で必要な用語が入っていないため，各用語を説明しながら当てはめていきます．また，今回の記事のメインである確率の必要性と位置付けについてもお話しします．

準備をまとめて実行しておきます．
```python:ダウンロードとインストール
!git clone https://github.com/ghmagazine/python_stat_sample.git
!mv ./python_stat_sample/data ./data
!pip install japanize-matplotlib
```
```python:インポートと設定
import numpy as np
import pandas as pd

import matplotlib.pyplot as plt
import japanize_matplotlib
import seaborn as sns

%precision 3
pd.set_option('precision', 3)
```

実行結果は省略．

### 母集団（Population）
観測したい対象データ全体のことを**母集団**と言います．今回のサンプルデータでは`400`人分の点数が母集団になります．Aさんにとって手に入れることができないデータですが，Aさんによる抽出をシミュレートするため，読み込みます．

```python:母集団の読み込み
all_scores = pd.read_csv('data/ch4_scores400.csv')['点数']
all_scores
```
```:実行結果
0      76
1      55
2      80
3      80
4      74
       ..
395    77
396    60
397    83
398    80
399    85
Name: 点数, Length: 400, dtype: int64
```

### 標本（Sample）
推測に使用する母集団の一部のデータのことを**標本**と言います．Aさんが手に入れた`20`人分の点数が，それに当たります．

### 標本抽出（Sampling）
母集団から標本を取り出すことを**標本抽出（Sampling）** と言います．Aさんが聞いて回って抽出したことが，それに当たります．

```python:標本抽出
sample_scores = np.random.choice(all_scores, 20, replace=False)
sample_scores
```
```:実行結果
array([55, 65, 69, 58, 63, 57, 86, 61, 50, 54, 64, 64, 65, 86, 48, 82, 75,
       50, 64, 58])
```
上記では`np.random.choice`によってランダムに標本を抽出する**無作為抽出（Random Sampling**）を行っています．また，`replace=False`としているので，同じ標本は一度しか抽出されない**非復元抽出**です．`replace=True`とすると同じ標本が再度抽出される可能性がある**復元抽出**になります．また，今回`20`人としている抽出数のことを**サンプルサイズ**と言います．

ここで，母集団から抽出したデータが偏ってしまう可能性についてに注意してください．例えば，Aさんの知り合い`20`人に点数を聞いた場合，性別や学年による偏りが出てしまうかもしれません．また，データが偏らないように校内で偶然見かけた人に声をかけたとしても，点数が良くなかった生徒は点数を教えてくれないかも，，，更に，プログラムで行ったように完全にランダムに抽出を行った場合でも，偶然データが偏ってしまうことが考えられます．

データの偏りは，抽出方法を工夫することで軽減できますが，完全には排除できません．推測統計では，このようにサンプルサイズや抽出方法で偏ってしまったデータを適切に扱うことが求められます！

### 標本統計量（Sampling statistic）
抽出した標本から計算される平均や分散などの統計量です．また，それぞれの統計量のことを標本平均，標本分散などと呼びます．

```python:平均と分散を求める
sample_scores.mean(), sample_scores.var()
```
```:実行結果
(63.700, 119.910)
```

### 母数（Parameter）
母集団の統計量（平均や分散など）のことを**母数**と言います．また，それぞれの統計量のことを母平均，母分散などと呼びます．母集団の全データを手元に得ることは難しいです．従って，母数は直接算出できないことが多いです．

### 推定量（Estimator）
例えば標本平均から母平均を推測するとき，標本平均は母平均の**推定量**である，と言います．また，推測された実際の値を**推定値**と言います．今回の例では，`400`人の母集団から`20`人を抽出し，抽出した標本から推定量として標本平均を計算して`63.7`という推定値が得られたことになります．

### 確率（Probability）
ここで，無作為抽出と統計量計算を繰り返してみます．具体的には，`400`人の点数から`20`人分の点数を無作為に選び，標本統計量を算出して表示する処理を`10`回繰り返しています．

```python:繰り返し標本抽出を実行
for i in range(10):
    sample_scores = np.random.choice(all_scores, 20, replace=False)

    mean, std = sample_scores.mean(), sample_scores.std()
    print(f'{i + 1:02}回  平均: {mean:.1f}  標準偏差: {std:.2f}')
```
```:実行結果
01回  平均: 66.7  標準偏差: 14.06
02回  平均: 72.9  標準偏差: 14.74
03回  平均: 76.3  標準偏差: 14.21
04回  平均: 71.2  標準偏差: 17.88
05回  平均: 73.6  標準偏差: 13.00
06回  平均: 64.4  標準偏差: 13.05
07回  平均: 72.5  標準偏差: 14.74
08回  平均: 73.0  標準偏差: 15.70
09回  平均: 70.8  標準偏差: 10.60
10回  平均: 68.7  標準偏差: 13.74
```
`20`人の標本抽出だと統計量が結構バラつきますね..標本抽出のところで述べたデータの偏りが影響しています．また，Aさんの例では`20`人の標本抽出を`10`回行うのは大変だと思われます..我々は限られた標本から得た推定値の曖昧さ（不確定さ）を定量的に認識する必要があります．

この不確定な現象を考えるときに有効なのが，**確率**です．確率とは偶然起こる現象の頻度ですが，一連の現象をモデル化して考えることで，今回の例を含む様々な現象の不確定さを定量的に表したり，制御できるようになります！

先程示した図を，ここで挙げた各用語に置き換えた全体像が以下です．

![スクリーンショット 2022-07-31 15.35.29.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/44273/2c0f86ce-f801-545f-7fed-cfe28b49c190.png)

図を文章にまとめてみました．
- 母集団全体を観測することができないため，母数は直接計算できない
- 母集団の一部を抽出した標本から標本統計量を計算し，母数を推定することができる
- 標本が偏ることで，推定には不確定さが伴う
- 確率を活用することで，不確定さを定量的に扱い制御することができる

## 確率モデルの全体像
先程挙げたような不確定な現象を数学的にモデル化したものを**確率モデル（Probability model**）と言います．これからサイコロを例に用いて確率モデルの全体像を説明していきます．モデル化については以下を参照．

https://kotobank.jp/word/%E6%95%B0%E5%AD%A6%E7%9A%84%E3%83%A2%E3%83%87%E3%83%AB-83261

ここで突然サイコロが例として出てきた理由は，確率モデルを考える上で都合が良く，シンプルで誰もが知っているツールだからです．また，世の中の様々な現象を単純なサイコロに置き換えて考えることができます！説明の中で随時例示していきます..

### 確率変数（Random Variable）
一般的なサイコロを投げたとき，出目は必ず1〜6の数字になります．一般的なサイコロは正方形なので，それぞれの数字が出る確率は同じになります．このように，試してみないと結果はわからないが，とりうる値と値が出る確率が決まっているものを，**確率変数**と言います．

また変数が出てきましたね..言葉だと理解し辛いですが，プログラムで定義される変数を想像すると分かりやすいかもしれません．以下はサイコロの出目が入る確率変数`p`を定義しています．
```python:1〜6が入る変数
p = np.random.randint(1, 7)
p
```
```:実行結果
2
```
サイコロの場合，1〜6がそれぞれ同じ確率で入る変数になります．サイコロ以外でも様々な確率変数を考えることができます．
| 説明 | 確率変数名 | とりうる値 |
|:--|:--|:--|
| サイコロ | 出目 |  `1`，`2`，`3`，`4`，`5`，`6` |
| コイン投げ | 向き |  `表`，`裏` |
| タイタニック | 生き残ったか | `生存`，`死亡`  |
| アプリレビュー | 星の数 | `1`，`2`，`3`，`4`，`5` |
| テスト | 点数 |  `0`〜`100` |
| ギャンブル | 残金 | `0`〜`∞` |
| ECサイト | 成果 | `離脱`，`カートに入れる`,`購入` |

ここで，確率が決まっていない現象は確率変数なのか，という疑問が浮かびます．例えばコイン投げであれば確率はほぼ半々だと分かっていますが，テストの点数はやってみないと結果が出ません．しかしよく考えてみると，実際にはコインの表面にも微妙に凹凸があって，確率を正確に算出するには無限回コインを投げる必要があります．逆にこれらの現象はすべて，無限回試せば確率を算出することができそうです．

つまり，無限回試したときに確率が確定する現象であれば，それを確率変数として考えることができるはず，，，と，私は解釈したのですが，ちょっと合っているか自信が無いです．正確な確率変数の定義は以下を参照ください．

https://ja.wikipedia.org/wiki/%E7%A2%BA%E7%8E%87%E5%A4%89%E6%95%B0

### 試行（Trial）
確率変数の結果を得ることが`試行`です．サイコロを投げると結果（出目）が分かるので，サイコロを投げることが試行に当たります．サイコロを10回投げた結果を得るプログラムが以下です．
```python:10回試行
[np.random.randint(1, 7) for _ in range(10)]
```

### 実現値（Realization）
試行によって観測された値です．サイコロを投げた結果，出目という実現値を得ます．つまり，↑の実行結果が実現値になります．

```:実行結果
[1, 1, 6, 5, 2, 3, 2, 4, 5, 2]
```

### 事象（Event）
試行の結果として起こりうる出来事を**事象**と言います．事象の例は以下．
- サイコロの出目が`6`
- サイコロの出目が`偶数`
- テストの点数が`92`点
- 大学受験に`合格`！
- Webサイトにアクセスがあったが`購入されず`..

出目が`6`は`1/6`，出目が`偶数`は`1/2`のように，確率はそれぞれの事象に対して対応づきます．また，これ以上分解できない出目`6`や`92`点などの事象を**根本事象（Elementary event**）と言います．更に，同時に起こり得ない事象を**排反事象（Exclusive event**）と言います．ここで重要な確率の性質として，**互いに排反**な事象において，そのうち少なくとも１つが起こる確率は，各事象の確率の和になる，というものがあります．いつだったか学校で習った**和の法則**ですね，，，

例えば，サイコロの出目`2`と出目`6`は同時に起こり得ないので互いに排反，と言えます．このときそれぞれの出目の確率は`1/6`なので，出目が`2`または`6`となる確率は，
```math
\frac{1}{6} + \frac{1}{6} = \frac{1}{3}
```
のように計算できます．逆に出目`2`と　出目`偶数`は，同時に起こり得るので排反ではなく，上記のような加算は成り立ちません．

　確率における和の法則や積の法則については以下

https://univ-juken.com/sekinohosoku-wanohosoku

### 確率分布（Probability distribution）
確率変数がどのような振る舞いをするのか表したものを**確率分布**と言います．どのような値がどのような確率で出現するか，とりうる値の分布を表します．まず，サイコロの確率分布を作ってみます．

```python:サイコロの確率分布
dice_df = pd.DataFrame(index=[1, 2, 3, 4, 5, 6])
dice_df['確率'] = 1 / 6
dice_df.T
```
||1|2|3|4|5|6|
|---|--:|--:|--:|--:|--:|--:|
|確率|0\.167|0\.167|0\.167|0\.167|0\.167|0\.167|

Markdownのフォーマットを少し修正しています..`0.167`というのは，約`1/6`です．確率分布を合計すると`1.0`になります．
```python:確率分布の合計
dice_df['確率'].sum()
```
```:実行結果
1.000
```
コイン投げについても同様に確率分布を作りました．
```python:コインの確率分布
coin_df = pd.DataFrame(index=['表', '裏'])
coin_df['確率'] = 1 / 2
coin_df.T
```
||表|裏|
|---|--:|--:|
|確率|0\.5|0\.5|

ここでサイコロとコインの違いについて考えてみましょう．

サイコロは根本事象が`6`種類，コインは`2`種類で，それぞれの根本事象の確率はサイコロが`1/6`で，コインが`1/2`です．他にも形や投げ方などに違いはありますが，確率モデルとして単純化して考えると，根本事象の数とそれに対応する確率が異なるだけだと分かります．

またサイコロには，コインと同じように`1`と`2`しか出目が無いもの，`12`の出目があるもの（正十二面体）など様々な種類があります．サイコロの種類については以下．思った以上に多種ですね，，，

https://oreore.red/polyhedron-dice-1-100

つまり，コインもサイコロの1つだと捉えることができます．更に，これまで例に挙げた世の中の物事について，何らかの形をしたサイコロを投げた結果だと考えることもできそうです．言い換えると，これまで例に挙げた様々な物事を確率モデルとして捉えることができる．そこから確率分布を導き出すことが，様々な物事を理解するのに役立ちそうです！

ただし，サイコロやコインのように最初から確率分布が判明しているケースは稀です．完全な確率分布を得るには，試行を無限回繰り返してそれぞれの根本事象に対応する確率を確定する必要がありますが，それは現実的ではありません..

### 度数分布と確率分布
そこで，試行を有限回繰り返して得られたデータから相対度数分布を作り，そこから確率分布を推測する，というアプローチを取ります．

相対度数については前記事の[こちら](https://qiita.com/o93/items/91feef63ea70070462d9#%E7%9B%B8%E5%AF%BE%E5%BA%A6%E6%95%B0relative-frequency)を参照ください．

このアプローチが有効であるかを確認するため，いかさまサイコロを用いた検証を行います．いかさまサイコロは重心がズレており，各出目（`1`〜`6`）がどの確率で出るかが偏っています．このいかさまサイコロを繰り返し投げて，確率分布を推測していきます！

```python:我々は知らない，いかさまサイコロの確率分布
roll_index = pd.Index([1, 2, 3, 4, 5, 6], name='出目')
loaded_dice_df = pd.DataFrame(index=roll_index)
loaded_dice_df['確率'] = [1/21, 2/21, 3/21, 4/21, 5/21, 6/21]
loaded_dice_df
```
|出目|確率|
|:-:|--:|
|1|0\.048|
|2|0\.095|
|3|0\.143|
|4|0\.190|
|5|0\.238|
|6|0\.286|

上記がいかさまサイコロの確率分布です．出目が大きいほど確率を高く設定しています．この確率分布は試行のシミュレーションに使うだけで，**我々にこの確率分布を知ることはできない**，とします．

いかさまサイコロ（確率変数）の試行には`np.random.choice`を使用し，上記の確率分布に従って試行を繰り返します．
```python:いかさまサイコロを100回試行
num_trial = 100
sample = np.random.choice(loaded_dice_df.index, num_trial, p=loaded_dice_df['確率'])
sample
```
```:実行結果
array([5, 5, 4, 2, 6, 5, 4, 2, 2, 6, 4, 6, 5, 6, 3, 6, 4, 6, 2, 4, 4, 3,
       5, 5, 4, 1, 5, 6, 4, 6, 4, 2, 3, 4, 4, 6, 6, 6, 4, 3, 6, 2, 2, 5,
       3, 1, 6, 5, 5, 5, 5, 6, 1, 6, 6, 5, 4, 3, 5, 6, 6, 6, 3, 5, 2, 5,
       6, 6, 5, 5, 4, 3, 2, 4, 4, 2, 5, 4, 6, 2, 3, 2, 3, 3, 5, 4, 4, 5,
       2, 5, 5, 4, 6, 4, 5, 5, 4, 3, 3, 6])
```
なんとなく出目が大きい数に偏っていてイカサマ臭が..

得られたサンプルデータから度数分布表を作ります．
```python:100回試行時の度数分布表
freq, _ = np.histogram(sample, bins=6, range=(1, 7))
pd.DataFrame({'度数': freq, '相対度数': freq / num_trial}, index=roll_index)
```
|出目|度数|相対度数|
|:-:|--:|--:|
|1|6|0\.06|
|2|9|0\.09|
|3|13|0\.13|
|4|14|0\.14|
|5|24|0\.24|
|6|34|0\.34|

`100`回試行したときの，それぞれの出目の数と割合が分かりました．いかさまサイコロの実際の確率分布と共に，`100`回試行したときの相対度数分布をヒストグラムにしてみます．

```python:100回試行時のヒストグラム
def show_hist_loaded_dice():
    plt.figure(figsize=(10, 6))
    plt.hist(sample, bins=6, range=(1, 7), density=True, label=f'{num_trial}回試行による相対度数分布')
    plt.hlines(loaded_dice_df['確率'], np.arange(1, 7), np.arange(2, 8), color='red', lw=3, label='実際の確率分布')
    plt.xticks(np.linspace(1.5, 6.5, 6), labels=np.arange(1, 7))
    plt.xlabel('出目')
    plt.ylabel('相対度数')
    plt.legend()
    plt.grid(True)
    plt.show()

show_hist_loaded_dice()
```
![スクリーンショット 2022-08-06 10.42.14.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/44273/2369423c-c5fb-f9e2-cfa8-daa2b0818a26.png)

実際の確率分布と相対度数分布を見比べると，傾向は似ているもののまだ差があるように思えます．このグラフの赤い線は，我々が見つけたい実際の確率分布であって，我々には赤い線は不明な状況です．実際の確率分布と，それを推測するために求めた相対度数分布との差を確認するため，敢えて正解データを表示していることに注意してください！

更に試行回数を`5000`件に増やしてみます．
```python:5000回試行時のヒストグラム
num_trial = 5000
sample = np.random.choice(loaded_dice_df.index, num_trial, p=loaded_dice_df['確率'])

show_hist_loaded_dice()
```
![スクリーンショット 2022-08-06 11.02.44.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/44273/a644a871-461a-2a1b-dac9-78af3d40d65c.png)
実際の確率分布に近くなってきましたね．この検証から，試行を繰り返して相対度数を求めることで，確率変数の確率分布を推測できそうだ，ということが分かりました．このまま試行回数を増やしていけば，相対度数分布は確率分布に収束していきます．

### 全体図
一通り確率モデルについて説明が終わりました．ここで確率モデルの全体図をまとめておきます．まずいかさまサイコロによるこれまでのアプローチを図示してみます．

![スクリーンショット 2022-08-06 12.15.51.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/44273/d12d3d57-2581-5f57-795f-2aac2c3a43cc.png)

いかさまサイコロを確率変数としたとき，その確率分布を直接知るには無限回の試行が必要で，実施は不可能です..そこで，サイコロを`100`回投げて，出目の記録を取って相対度数分布を求めました．そこからいかさまサイコロの真の確率分布を推定しよう，というアプローチです．

いかさまサイコロを一般的な用語に置き換えたものが以下です．

![スクリーンショット 2022-08-07 7.57.56.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/44273/f2690bea-6dfd-873f-8fd0-eb3be545f17b.png)

アプリレビューやWebサイトへのアクセスなど，様々な現象を確率変数として考えることで，同様のアプローチが可能になります．また，この図と推測統計の図は構成が似ています，，，どちらも左側の要素のすべてを観測することはできない，ということが大前提になっています．今回参考にしている書籍を読んだタイミングではないのですが，ここからの個人的な学びとして，

- **事実**は**真実**ではない
- **事実**を**解釈**して**真実**を**推定**するしかない

という思想を得ました．上図において，事実は右側，真実は左側です．これは拡大解釈し過ぎかもしれませんが，機械学習やデータ分析だけでなく，自分の人生にも影響を及ぼしていたりします..

## 推測統計と確率モデル
基本の最後に，推測統計の説明に用いた点数の無作為抽出の例を確率モデルで考えてみます．母集団から無作為抽出によって得られる標本を，母集団の確率分布に従う確率変数と見なします．ただし，母集団のすべてが観測できない（無限回試行できない）限り，正確な確率分布は算出できません．

まずは全生徒`400`人分の点数の相対度数分布をヒストグラムで表示します．
```python:全生徒のヒストグラム（相対度数分布）
def show_hist_scores(s, bins=100):
    plt.figure(figsize=(10, 6))
    plt.hist(s, bins=bins, range=(0, 100), density=True)
    plt.xlim(20, 100)
    plt.xlabel('点数')
    plt.ylabel('相対度数')
    plt.grid(True)
    plt.show()

all_scores = pd.read_csv('data/ch4_scores400.csv')['点数']
show_hist_scores(all_scores)
```
![スクリーンショット 2022-08-07 9.01.54.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/44273/f876caa1-7ed2-ef15-7d12-fa3646d850fa.png)

この度数分布を母集団の確率分布だと見なすと，例えば`70`点の人は全体の`2.5`%を占めているので，無作為抽出を`1`回試行したとき，`70`点を得る確率は`2.5`%だということが分かります．次に試行をシミュレートしてみます．

```python:無作為抽出の試行1回
np.random.choice(all_scores)
```
```:実行結果
76
```
`76`点という実現値を得ました．

では，`100`回試行を繰り返して，同様にヒストグラムを表示してみます．
```python:試行100回のヒストグラム
scores = np.random.choice(all_scores, 100)
show_hist_scores(scores)
```
![スクリーンショット 2022-08-07 12.24.40.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/44273/78b9f779-be81-1c72-3e1e-a58f9207d0f8.png)

荒いですが，同じ傾向は見受けられます．今度は試行を`10000`回にしてみます．
```python:試行10000回のヒストグラム
scores = np.random.choice(all_scores, 10000)
show_hist_scores(scores)
```
![スクリーンショット 2022-08-07 12.33.13.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/44273/9f520fb2-fe6c-d0db-b923-3287ac570432.png)
母集団の確率分布に近づいてきました，，，

ここで注意したいのは，今回のケースでは確率変数が`400`人のテストの点数だという点です．仮にこのテストを日本中の高校生が受けたとしたら，確率分布は上記とは異なるものになるはず，，，我々が対象とする確率変数や母集団をどう定義するかで，推定すべき確率分布が変わる，ということです．

推測統計と確率モデルの関係をまとめます．
- 推測統計は確率モデルを使って母集団の統計的性質を推測する手法
- 母集団から得られた標本1つ1つが，確率分布に従う確率変数だと見なせる
- 母集団からの無作為抽出は，確率分布に従う確率変数の試行に当たる
- データは結果が確定した値なので，実測値に該当する

母集団を既に結果が確定しているデータと捉えるのではなく，確率分布に従って毎回結果を生み出す**歪な形のサイコロ**だとイメージすると，理解しやすいかもしれません．抽象化され過ぎていて「物事はそんな単純じゃない！」という批判が聞こえてきそうですが，，，物事には色んな捉え方があって，それが有用なら使えば良いんじゃないかと，個人的には思っています．

## 統計量の確率変数
最後に，統計量を確率変数として捉える例を示します．

**Aさんが無作為に`20`人に声をかけて得た点数の平均を算出**

これを`1`回の試行として捉えると，`20`人の標本を抽出して算出した平均値を，確率変数の実現値だと考えることができます．それでは，この確率変数が従う確率分布を推定してみましょう．`20`人分の点数の抽出を`10000`回シミュレートして相対度数分布を表示しています．

```python:20人分の標本平均の相対度数分布
sample_means = [np.random.choice(all_scores, 20).mean() for _ in range(10000)]

plt.figure(figsize=(10, 6))
plt.hist(sample_means, bins=100, range=(0, 100), density=True, label='相対度数分布')
plt.vlines(np.mean(all_scores), 0, 1, color='red', label='母平均')
plt.xlim(50, 90)
plt.ylim(0, 0.15)
plt.xlabel('点数')
plt.ylabel('相対度数')
plt.legend()
plt.grid(True)
plt.show()
```
![スクリーンショット 2022-08-13 9.14.39.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/44273/cc72078b-50b9-c6b2-74a0-a5b83ca43305.png)
ばらついているものの，母平均を中心に求めた平均値が分布していることから，このシミュレートによって標本平均から母平均が推定できることが分かります．

また，`20`人の標本抽出だとこのくらい平均値がばらつく，ということが読み取れます．このばらつきは，今実施している分析手法の曖昧さを定量的に表しています．確率モデルを上手く使えば，母集団の統計的性質を推定するだけでなく，分析によって算出された値の曖昧さも定量的に捉えることができます！

# 離散型確率変数（Discrete random variable）
ここまでは推測統計と確率の基本に触れてきました．ここからは，確率変数の種類や性質についてより詳細に学んでいきます．確率変数はいくつかの種類に分類できます．最初は**離散型確率変数**から．サイコロの出目やレビューの★の数のように，とりうる値が離散している確率変数です．

ちなみに書籍ではこの辺りから数式の記載が増えています．記事ではなるべく数式は記載せず，プログラムと文章での説明を行っています．繰り返しになりますが，数式が気になる方は書籍を購入ください！

## 1次元
最初に1次元の確率変数から説明していきます．次元と言われると難しく感じるかもしれませんが，1つの確率変数がある，と考えればOKです．例えば1つのサイコロを，1次元の離散型確率変数，と捉えることができます．

### 確率質量関数／確率関数（Probability mass function／PMF）
1次元の離散型確率変数を定義するために必要なのは，確率変数がとりうる値と，その値が出現する確率です．2つの値の対応であればPythonの`dict`を使えば解決しそうですが，今後この対応が複雑になることを考えて，関数で定義することにします．これを**確率質量関数**，または**確率関数**と言います．とりうる値の1つを引数として入力したら，確率分布に従った確率を出力する関数です．

いかさまサイコロを例に，最初にとりうる値を作ります．サイコロの出目です．
```python:確率変数がとりうる値
x_set = np.array([1, 2, 3, 4, 5, 6])
```
次に確率関数を定義します．
```python:確率関数
def f(x):
    return x / 21 if x in x_set else 0
```
引数`x`が`x_set`内にあるとき， `x / 21`を確率として返しています．パフォーマンス悪くなりそうですが，`21`を`np.sum(x_set)`としても良さそう．

確率変数は，とりうる値と確率のセットなので，`X`として作っておきます
```python:確率変数
X = [x_set, f]
```
大文字で変数に名前を付けるのは違和感がありますが，機械学習のコードだと良く見かけますね．とりあえずその辺りは書籍に合わせています..

確率関数を使用して確率分布を作ります．
```python:確率分布
prob = np.array([f(x) for x in x_set])
dict(zip(x_set, prob))
```
```:実行結果
{1: 0.048, 2: 0.095, 3: 0.143, 4: 0.190, 5: 0.238, 6: 0.286}
```
出目が大きくなればなるほど出現確率が増える，いかさまサイコロの確率分布を作ることができました！

そして確率は，すべて`0.0`以上で，総和は`1.0`になる性質を持っています．
```python: 確率の性質
print(f'すべて0以上: {np.all(prob >= 0)}')
print(f'総和: {np.sum(prob):.3f}')
```
```:実行結果
すべて0以上: True
総和: 1.000
```
この辺りから，これらが何の役に立つのか解らなくなってきますが，，，例えば営業の電話をかける「架電」を確率変数とした場合，
- 出る
- 出ない
- 話し中
- 留守電
- 未使用の番号

というような，とりうる値（事象）が考えられます．様々な時間帯や曜日，季節における架電結果（確率分布）を知ることは，効率良く架電を行う上で重要な情報になります．試行をシミュレートするための確率変数を定義したり，確率の性質を知ることは，できるだけ少ない試行で架電の確率分布を推定するのに役立つはず..

難しい用語や定義からは用途をイメージし辛いですが，そんなモチベーションで読み進めて頂ければと思います！

### 累積分布関数／分布関数（Cumulative distribution function／CDF）
確率変数`X`が引数`x`以下になる確率です．とりうる値が`x`以下の確率を加算して返す関数です．

```python:累積分布関数
def F(x):
    return np.sum([f(xx) for xx in x_set if xx  <= x])

F(4)
```
```:実行結果
0.476
```
いかさまサイコロを投げたとき`4`以下が出る確率は，`47.6`%だ，ということです．

2つの関数から得られる分布は，前記事で挙げた，有限のデータから得られる相対度数分布，累積相対度数分布に似ています．相対度数分布は試行を繰り返すことで確率分布に収束していく，同様に累積相対度数分布は累積の確率分布に収束していくイメージです．

### 確率変数の変換
ここで，確率変数のとりうる値を置き換える操作について説明します．いかさまサイコロの出目に対して`10X - 10`の報酬を与える，とした場合，確率変数`X`を変換して別の確率変数`Y`を作ることができます．例えば出目が`1`なら報酬は`0`，出目が`5`なら`40`になります．この確率変数の確率分布を作ってみます．

```python:とりうる値を置き換えて確率分布を作成
y_set = np.array([10 * xx - 10 for xx in x_set])
dict(zip(y_set, prob))
```
```:実行結果
{0: 0.048, 10: 0.095, 20: 0.143, 30: 0.190, 40: 0.238, 50: 0.286}
```
`x_set`から`y_set`を作って置き換えた上で，いかさまサイコロの確率は変わらないため，`prob`はそのまま流用しています．

報酬を変更した上で次項で出てくる期待値を求めるときや，確率変数の標準化（平均を引いて標準偏差で割る）に必要な操作らしいのですが，現状の情報だけではいまいちピンときませんね..少し調べてみたのですが，この段階では，確率変数から別の確率変数を作ったり，確率変数同士の演算ができるよ，くらいの理解で良さそうです．

確率変数の標準化については以下

https://k-san.link/standardized/

### 期待値（Expected value）
確率変数を無限回試行して得られた実現値の平均を**期待値**と言います．無限回試行することはできませんが，確率変数のとりうる値と確率の積の総和で算出することができます．有限回の繰り返しで得られたデータの平均は，試行を増やしていくと期待値に収束する，ということです．つまり，標本から算出される平均・分散などの指標は，標本が増加するにつれて確率変数の指標に収束していきます！

では，いかさまサイコロの確率変数`X`の期待値を計算してみます．

```python:いかさまサイコロXの期待値
def E(X, g=lambda x: x):
    x_set, f = X
    return np.sum([g(xx) * f(xx) for xx in x_set])

E(X)
```
```:実行結果
4.333
```
引数に確率変数`X`を入力して，とりうる値と確率の積を算出，`np.sum`で総和しています．引数には`g`という関数を指定できます．デフォルト値はとりうる値をそのまま使用しますが，先程説明した確率変数を変換する関数などを指定することができます．先程の確率変数`Y　=　10X - 10`を指定してみます．

```python:確率変数Yの期待値
E(X, g=lambda x: 10 * x - 10)
```
```:実行結果
33.333
```
いかさまサイコロの出目ではなく，出目から計算される報酬の期待値は`33.333`になりました．

また，期待値は**線形性**が成り立ちます．

```python:線形性の確認
10 * E(X) - 10
```
```:実行結果
33.333
```
`E(X)`を求めてから変換を行っても同じ結果になりましたね！関数の線形性については以下

https://cognicull.com/ja/g4y2mgcl

### 分散（Variance）
データにおける分散と同様に確率変数のばらつきを表します．前記事の分散の記事は[こちら](https://qiita.com/o93/items/91feef63ea70070462d9#%E5%88%86%E6%95%A3variance)．では，いかさまサイコロの確率変数`X`の分散を算出してみます．

```python:いかさまサイコロXの分散
def V(X, g=lambda x: x):
    x_set, f = X
    mean = E(X, g)
    return np.sum([(g(xx) - mean) ** 2 * f(xx) for xx in x_set])

V(X)
```
```:実行結果
2.222
```
`g(xx) - mean`が偏差に当たります．それを`2`乗し，とりうる値の確率`f(xx)`との積を算出，更に`np.sum`で総和し，確率変数の分散を求めています．期待値のときと同様に，確率変数Yの分散`Y　=　10X - 10`を指定してみます．

```python:確率変数Yの分散
V(X, g=lambda x: 10 * x - 10)
```
```:実行結果
222.222
```

また，期待値のときは線形性が成り立ちましたが，分散の場合は**分散の公式**が成り立つとのこと．

```python:分散の公式
10 ** 2 * V(X)
```
```:実行結果
222.222
```
係数`10`（傾き）を`2`乗するだけで，切片`-10`は式から消えていますが，結果は一致しました！公式の説明は割愛しますので，以下を参照ください．

https://risalc.info/src/st-expectation-variance-properties.html

## 2次元
続いて２次元の離散型確率変数について説明します．2次元というのは，2つの確率変数を同時に扱うということです．前記事の[2次元データ](https://qiita.com/o93/items/91feef63ea70070462d9#2%E6%AC%A1%E5%85%83%E3%83%87%E3%83%BC%E3%82%BF)に似ていますね．1次元のときと同様に，データが無限に存在する（無限回の試行）と考えたときの共分散などが，2つの確率分布から算出できそうです！

### 同時確率分布／同時分布（Joint probability distribution）
２つの確率変数の振る舞いを同時に考えた分布のことを**同時確率分布**と言います．2つのいかさまサイコロを同時に投げても，一方の出目がもう一方の出目に影響しません．それでは共分散が`0.0`になってしまって都合が悪いため，サイコロAの出目を`Y`，サイコロAの出目とサイコロBの出目を足したものを`X`として，わざと2つの変数に関連を持たせます．

まとめると，
- サイコロAとサイコロBを同時に投げて出目を確認
- `X`・・・サイコロAの出目＋サイコロBの出目
- `Y`・・・サイコロAの出目

を試行とします．まず，2次元の確率変数がとりうる値を考えます．
```python:とりうる値
x_set = np.arange(2, 13)
y_set = np.arange(1, 7)

print(f'Xのとりうる値: {x_set}')
print(f'Yのとりうる値: {y_set}')
```
```:実行結果
Xのとりうる値: [ 2  3  4  5  6  7  8  9 10 11 12]
Yのとりうる値: [1 2 3 4 5 6]
```
`x_set`は2つのサイコロの出目の和なので`2〜12`，`y_set`は`1〜6`になります．この2つの配列の組み合わせが，今回の確率変数のとりうる値です．ではそれぞれの確率はどうなるでしょうか，，，いかさまサイコロの確率分布は以下でした．

|出目|確率|
|:-:|--:|
|1|1 / 21|
|2|2 / 21|
|3|3 / 21|
|4|4 / 21|
|5|5 / 21|
|6|6 / 21|

この確率分布に従うサイコロA，Bを同時に投げたとき，Aの出目が`Y`，AとBの和が`X`なので，`X`からサイコロAの出目を引けば，`X`のサイコロBの出目が分かります．あとはサイコロAの出目の確率とサイコロBの出目の確率の積が，`X`と`Y`を組み合わせた確率になります．学校で習った積の法則ですね！

例えば，`X`が`5`，`Y`が`2`のとき，サイコロAの出目は`2`，サイコロBの出目は`5 - 2 = 3`なので，確率は
```math
\frac{2}{21} \times \frac{3}{21} = \frac{6}{441}
```
と計算できます．サイコロA，Bの出目から`X`と`Y`を求めて，そこから確率を算出するほうが直感的ですが，今回サイコロA，Bの出目を加算するところは試行の中で行われていることなので，あくまで試行結果からサイコロBの出目を逆算することにしています．

では，確率分布を作ってヒートマップを表示してみます．

```python:確率分布を表示
def f_XY(x, y):
    a, b = y, x - y
    if 1 <= a <= 6 and 1 <= b <= 6:
        return  a * b / (21 * 21)
    else:
        return 0

prob = np.array([[f_XY(i, j) for j in y_set] for i in x_set])
prob_df = pd.DataFrame(prob, index=pd.Index(x_set, name='X'), columns=pd.Index(y_set, name='Y'))

plt.figure(figsize=(8, 7))
sns.heatmap(prob_df, annot=True, fmt='1.3f', cmap='Reds')
plt.show()
```
![スクリーンショット 2022-08-14 7.38.19.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/44273/25ef4b45-7db2-1f72-caec-42affef501de.png)
サイコロは`1〜6`までしか出目が無いため，例えば$Y$が`1`のとき$X$は`8`以上になりません．そういったあり得ない確率はすべて`0.0`になっていますね．そして$X$，$Y$両方の数字が大きくなるにつれて確率が高くなっています．いかさまサイコロらしい傾向です．

2次元でも確率の性質を確認しておきます．
```python: 確率の性質
print(f'すべて0以上: {np.all(prob >= 0)}')
print(f'総和: {np.sum(prob):.3f}')
```
```:実行結果
すべて0以上: True
総和: 1.000
```
1次元のときと同様の結果が得られましたね．

### 同時確率関数（Joint probability function）
2つのとりうる値を入力として，その確率分布に従った確率を返す関数です．先程定義した`f_XY`ですね．改めて記載して説明します．

```python:同時確率関数
def f_XY(x, y):
    a, b = y, x - y
    if 1 <= a <= 6 and 1 <= b <= 6:
        return  a * b / (21 * 21)
    else:
        return 0
```
引数`x`，`y`からサイコロA，Bの出目`a`，`b`を求めて，それぞれがサイコロの出目の範囲に収まっているとき，積の法則に従って確率を算出しています．

確率変数`XY`も定義しておきます．
```python:確率変数XY
XY = [x_set, y_set, f_XY]
```

### 周辺確率分布／周辺分布（Marginal probability distribution）
複数の物事から成る確率分布に対して，1つの物事による確率分布を**周辺確率分布**と言います．今回の2次元の確率変数だと，`X`や`Y`それぞれの確率分布が周辺確率分布です．例えば2次元の確率変数`XY`から確率変数`X`を取り出すには，確率変数`Y`の影響を取り除く必要があります．では，`XY`から`X`と`Y`を取り出してみましょう．

```python:周辺分布Xの抽出
prob_df.sum(axis=1).to_frame().T
```
|2|3|4|5|6|7|8|9|10|11|12|
|--:|--:|--:|--:|--:|--:|--:|--:|--:|--:|--:|
|0\.002|0\.009|0\.023|0\.045|0\.079|0\.127|0\.160|0\.172|0\.166|0\.136|0\.082|

`X`のとりうる値それぞれで`Y`のとりうる値の確率を合計すれば，周辺分布を求めることができます．`DataFrame.sum`を使用しました．`Y`についても同様に求めます．

```python:周辺分布Yの抽出
prob_df.sum(axis=0).to_frame().T
```
|1|2|3|4|5|6|
|--:|--:|--:|--:|--:|--:|
|0\.048|0\.095|0\.143|0\.190|0\.238|0\.286|

上記と同じ処理を行う関数を用意して確率変数`X`と`Y`を作ります．
```python:それぞれの周辺分布に従う確率変数X，Y
def f_X(x):
    return np.sum([f_XY(x, y) for y in y_set])

def f_Y(y):
    return np.sum([f_XY(x, y) for x in x_set])

X = [x_set, f_X]
Y = [y_set, f_Y]
```
`f_XY`を使用して，引数ともう一方のすべてのとりうる値の確率を合計しています．作った関数を使用して，それぞれの周辺分布を表示してみます！

```python:周辺分布の表示
prob_x = np.array([f_X(k) for k in x_set])
prob_y = np.array([f_Y(k) for k in y_set])

fig = plt.figure(figsize=(16, 6))
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)

ax1.bar(x_set, prob_x)
ax1.set_title('Xの周辺分布')
ax1.set_xlabel('Xのとりうる値')
ax1.set_ylabel('確率')
ax1.grid()

ax2.bar(y_set, prob_y)
ax2.set_title('Yの周辺分布')
ax2.set_xlabel('Yのとりうる値')
ax2.set_ylabel('確率')
ax2.grid()

plt.show()
```
![スクリーンショット 2022-08-14 10.06.52.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/44273/f9f69ebd-b5f6-a5e1-41b0-1ed4621509f2.png)

`X`の周辺分布は2つのサイコロの足した数なので，`2`，`12`などは出目の組み合わせが少なく，確率が低くなります．しかし，使用しているのがいかさまサイコロなので，若干大きい数に確率の山が寄っているのが分かりますね，，，また`Y`はいかさまサイコロの確率分布そのものになっていることが分かります．

### 期待値（Expected value）
2次元の確率変数の期待値は1次元と同じです．関数`E2`を定義して`X`と`Y`それぞれの期待値を求めてみます．

```python:X，Yの期待値
def E2(XY, g):
    x_set, y_set, f_XY = XY
    return np.sum([[g(i, j) * f_XY(i, j) for i in x_set] for j in y_set])

mean_X = E2(XY, lambda x, y: x)
mean_Y = E2(XY, lambda x, y: y)
mean_X, mean_Y
```
```:実行結果
(8.667, 4.333)
```
関数`g`でとりうる値そのものを指定して返せば，期待値を求めることができます．例えばここにサイコロの出目に基づく報酬などを設定できるのは，1次元のときと同じですね．`x + y`とすれば`X`と`Y`の期待値の合計が算出できます，
```python:期待値の合計
E2(XY, lambda x, y: x + y)
```
```:実行結果
13.000
```
1次元のときと同様に期待値には線形性があるため，`X`と`Y`それぞれの期待値を`E2`で算出してから加算しても，結果は同じになります．
```python:期待値の線形性
mean_X + mean_Y
```
```:実行結果
13.000
```

### 分散（Variance）
分散も1次元と同様に，偏差の2乗の期待値によって求められます．

```python:X，Yの分散
def V2(XY, g):
    x_set, y_set, f_XY = XY
    mean = E2(XY, g)
    return np.sum([[(g(i, j) - mean) ** 2 * f_XY(i, j) for i in x_set] for j in y_set])

var_X = V2(XY, lambda x, y: x)
var_Y = V2(XY, lambda x, y: y)
var_X, var_Y
```
```:実行結果
(4.444, 2.222)
```

### 共分散（Covariance）
共分散は`X`と`Y`それぞれの偏差の積の期待値によって求められ．2つの確率変数にどの程度相関があるかを示します．[データ編の共分散](https://qiita.com/o93/items/91feef63ea70070462d9#%E5%85%B1%E5%88%86%E6%95%A3covariance)に似ていますね．

```python:X，Yの共分散
def cov(XY):
    x_set, y_set, f_XY = XY
    mean_X = E2(XY, lambda x, y: x)
    mean_Y = E2(XY, lambda x, y: y)
    return np.sum([[(i - mean_X) * (j - mean_Y) * f_XY(i, j) for i in x_set] for j in y_set])

cov_XY = cov(XY)
cov_XY
```
```:実行結果
2.222
```

### 相関係数（Correlation coefficient）
最後に相関係数です．こちらも[データ編の相関係数](https://qiita.com/o93/items/91feef63ea70070462d9#%E7%9B%B8%E9%96%A2%E4%BF%82%E6%95%B0correlation-coefficient)と同様で，共分散を`X`，`Y`の標準偏差で割って求めます．

```python:相関係数
def corr(XY):
    cov_xy = cov(XY)
    var_x = V2(XY, lambda x, y: x)
    var_y = V2(XY, lambda x, y: y)

    return cov_xy / np.sqrt(var_x * var_y)

corr(XY)
```
```:実行結果
0.707
```
`Y`はサイコロAの出目，`X`はサイコロAとBの出目を足したものなので，２つの確率変数には相関が仕組まれています．算出された相関係数は，`0.707`と，高い相関を示していることがわかります．

# おわりに
ここまでを確率編1とします．今回は，推測統計と確率モデルの全体像，1次元と2次元の離散型確率変数について記載しました．

我々の目的に必要な全データと，手元にあるデータが異なる場合，そのデータは母集団から抽出された標本であり，有限回の試行から生み出された実現値の集まりです．無限回の試行（全データ）から得られた確率分布を知ることができれば，我々は目的を果たすことができます．

目的に対するアプローチとして，手元にあるデータからなるべく正確に全データの性質を知ること，全データに対して手元にあるデータの曖昧さを定量的に知ることに，確率モデルが役立つ，ということを学びました．また，データ(有限の試行結果)から得られる統計値は，データを増やしていく（無限の試行結果に近づける）と，データを確率変数と捉えたときの統計値に収束していくことが分かりました．

個人的には，ふわっとした理解しか出来ていなかった確率モデルや，記述統計と推測統計の関係性について整理する良い機会になりましたが，この分野についての実務経験はあまり多くないので，間違っている点などあればご指摘頂ければと思います！最後まで読んで頂きありがとうございました！！

次回は代表的な確率分布や，連続型確率変数について記載しています．

https://qiita.com/o93/items/a2e02cf98126f47cd66c