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

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

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部構成の3つ目です．以下も参照ください．

- [準備〜データ編](https://qiita.com/o93/items/91feef63ea70070462d9)・・・投稿済み
- [確率編1](https://qiita.com/o93/items/1b3f379da3e62c6252d8)・・・投稿済み
- 確率編2・・・この記事
- 確率編3・・・まだ
- 推測編・・・まだ

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

今回もエンジニアに馴染みの無い用語が沢山出てきます．また，深堀りすればするほど，この技術が何に使えるのか？を見失ってしまい，読み進めるのを諦めてしまいがちです，，，そうならぬよう，書籍を参考にしつつも，補足を入れたり，情報の取捨て選択をしているので，書籍に沿って進めたい方は書籍を購入ください！

確率編2は，確率編1にて学習した離散型確率分布の詳細と，連続型確率分布について一通り学んでいきます．

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

準備をまとめて実行しておきます．今回新たに，科学計算ライブラリである[SciPy](https://scipy.org/)を使用します．

```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
from scipy import stats
from scipy.special import comb
from scipy.special import factorial
from scipy import integrate
from scipy import optimize
from functools import partial

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

import warnings
warnings.filterwarnings('ignore', category=integrate.IntegrationWarning)

%precision 3
pd.set_option('precision', 3)
pd.options.display.float_format = '{:.3f}'.format
```

# おさらい

限られた標本（データ）から母集団（ターゲット）の統計的特徴を推定し，目的を達成するために活用することが，我々が確率を学ぶモチベーションです．母集団から抽出できる標本が限られる場合，推定する統計的特徴は不確かです．そこで，**統計的特徴の推定**というタスクを，確率変数の振る舞いの推定，つまり**確率分布の推定**，というアプローチで進めることにしました．

確率変数がとりうる値を入れると，値の出現確率を返してくれる関数が**確率関数**です．例えば一般的なサイコロであれば，以下のように確率関数から確率分布を算出できます．

```python:サイコロの確率関数と確率分布
def dice_prob_func(x):
    return 1 / 6

dice_x_set = [1, 2, 3, 4, 5, 6]
dice_prob = np.array([dice_prob_func(x) for x in dice_x_set])

plt.bar(dice_x_set, dice_prob)
plt.ylim(0.0, 1.0)
plt.xlabel('とりうる値')
plt.ylabel('確率')
plt.grid(True)
plt.show()
```
![スクリーンショット 2022-09-10 11.17.06.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/44273/80be9048-2b00-591b-4934-80c30c05b4de.png)
一般的なサイコロは．出目の確率が`1/6`になることが明白です．しかし，世の中の様々な物事は確率分布が分かっておらず，これを推定することこそが我々の目的であるわけです．

```python:相対度数分布による確率分布の推定
num_trial = 100
dice_sample = np.random.choice(dice_x_set, num_trial)

plt.hist(
    dice_sample, bins=len(dice_x_set), range=(1, 7), density=True,
    label=f'{num_trial}回試行による相対度数分布')
plt.hlines(
    dice_prob, 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.ylim(0.0, 1.0)
plt.xlabel('とりうる値')
plt.ylabel('相対度数')
plt.legend()
plt.grid(True)
plt.show()
```
![スクリーンショット 2022-09-10 12.26.47.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/44273/8dbb6094-39e6-6620-6968-41bf376fbb33.png)
上記は，サイコロを`100`回投げたときの相対度数分布です．試行を増やすと，相対度数分布は確率分布に収束していきますが，沢山試行する（沢山のデータを得る）のは大変，，，少ない試行で効率良く確率分布を推定する必要があります．

# 代表的な離散型確率分布
先程挙げた，サイコロを投げた結果から確率分布を推定する方法は，予め確率分布の形状を仮定していません．このような手法を**ノンパラメトリック**な手法，といいます．この手法は，試行を沢山繰り返すことで確率分布を推定することができますが，あまり効率的ではありません..

それに対して，例えばサイコロの形状が正方形であることから，`1〜6`の出目の確率が大体均等になるだろう，と確率分布のおおよその形状を仮定してから試行を繰り返せば，より少ない試行で確率分布を推定することができます．このような手法を**パラメトリック**な手法，といいます．

具体的には，試行結果のみに基づいて確率関数を実装するのではなく，まず母集団の性質から確率分布の形状を仮定した上で確率関数を実装します．その上で，確率関数のパラメータを調整して，母集団の確率分布に近づけていきます．パラメータの調整については後日公開予定の記事（推測編）に記載する予定です．

まとめると，以下のようになります．
- ノンパラメトリックな手法
    - 母集団から確率分布の形状を仮定せず確率分布を推定
    - 形状を仮定しない分，確率分布を柔軟に表現できる
    - 正確な確率分布を得るには，沢山の試行結果が必要
- パラメトリックな手法
    - 母集団の性質から確率分布の形状を仮定し，少数のパラメータを調整して確率分布を推定
    - より少ない情報から確率分布を推定することができる
    - 形状を決め打ちするので，確率分布の表現力はノンパラメトリックに劣る
    - 仮定を誤ると確率分布の正確に推定できない可能性が..

書籍の大部分は，後者のパラメトリックな手法についての説明です．パラメトリックな手法を用いるには，**どのような場面でどのような確率分布を使用すべきか**を知る必要があります．ここからは，代表的な離散型確率分布と，それらがどのような場面で使われるかを学んでいきます！

## 準備

最初に離散型確率分布の性質を確認するための関数をいくつか用意します．

以下は前記事で学んだ期待値と分散を求める関数です．
```python:期待値と分散の関数
def E(X, g=lambda x: x):
    x_set, f = X
    return np.sum([g(k) * f(k) for k in x_set])

def V(X, g=lambda x: x):
    x_set, f = X
    mean = E(X, g)
    return np.sum([(g(k) - mean) ** 2 * f(k) for k in x_set])
```

以下は確率変数の性質を満たすか確認した上で，期待値と分散を表示する関数です．
```python:確率変数の性質の確認・表示
def check_prob(X):
    x_set, f = X

    prob = np.array([f(k) for k in x_set])
    assert np.all(prob >= 0), '負の確率があります'

    prob_sum = np.round(np.sum(prob), 6)
    assert prob_sum == 1, f'確率の和が{prob_sum}になりました'

    print(f'期待値: {E(X):.4}, 分散: {V(X):.4}')
```

最後は確率分布と期待値をグラフで表示するための関数です．
```python:確率分布・期待値・分散のグラフ表示
def plot_prob(X, size=(10, 6)):
    x_set, f = X

    prob = np.array([f(k) for k in x_set])
    e = E(X)

    plt.figure(figsize=size)
    plt.bar(x_set, prob, label='確率')
    plt.vlines(e, 0, 1, label='期待値', color='red')
    plt.xticks(np.append(x_set, e))
    plt.ylim(0, prob.max() * 1.2)
    plt.legend()
    plt.grid(True)
    plt.show()
```

## ベルヌーイ分布（Bernoulli distribution）
とりうる値が`0`と`1`しかない，最も基本的な離散型確率分布です．ベルヌーイ分布に従う確率分布の試行を**ベルヌーイ試行**と言い，実現値`1`を**成功**，実現値`0`を**失敗**と言います．とりうる値が2つしかないため，`1`が出る確率を`p`とした場合，もう片方の`0`が出る確率は`1 - p`と定まります．つまり，ベルヌーイ分布の形状を調整するために必要なパラメータは，成功（`1`が出る）確率である`p`だけです．

では早速，ベルヌーイ分布の確率関数を作ってみましょう．まずNumPyによる実装，次にSciPyをによる実装を確認していきます．

### NumPyによる実装

```python:ベルヌーイ分布の確率関数
def prob_bernoulli(x, p, x_set):
    if x in x_set:
        return (p ** x) * ((1 - p) ** (1 - x))
    else:
        return 0
```
引数は3つで，値`x`，先程挙げたパラメータ`p`，とりうる値`x_set`です．ベルヌーイ分布の場合，`x_set`は`[0, 1]`と決まっていますが，関数外でも使用するため，今回は敢えて引数にしています．`x_set`に`x`が含まれない場合は`else`側で`0`を返します．`if`側のコードは無駄にややこしく書いている印象を受けましたが..ベルヌーイ分布の定義をそのまま実装するとこうなるようです．プログラムの詳細を説明します．

例えば`p`が`0.2`のとき，プログラムは以下のように動作します．
- 値`x`が`1`のとき
    - `p ** x`・・・`0.2`の`1`乗は`0.2`
    - `(1 - p) ** (1 - x)`・・・`0.8`の`0`乗は`1`
    - `0.2 * 1 = 0.2`・・・`0.2`を返す
- 値`x`が`0`のとき
    - `p ** x`・・・`0.2`の`0`乗は`1`
    - `(1 - p) ** (1 - x)`・・・`0.8`の`1`乗は`0.8`
    - `1 * 0.8 = 0.8`・・・`0.8`を返す

`1`乗と`0`乗を上手く切り替えることで，出力を制御しています．`prob_bernoulli`を実行してみます．

```python:ベルヌーイ関数の実行
bern_x_set = np.array([0, 1])
bern_p = 0.2

bern_prob = np.array(
    [prob_bernoulli(k, bern_p, bern_x_set) for k in bern_x_set])
pd.DataFrame(
    bern_prob, columns=['確率'], index=pd.Index(bern_x_set, name='とりうる値'))
```
|とりうる値|確率|
|--:|--:|
|0|0\.8|
|1|0\.2|

1つ目が実現値`0`の確率で`0.8`，2つ目が実現値`1`の確率で`0.2`になっていますね！もう少し使い勝手を良くするため，ベルヌーイ分布に従う確率変数を返す関数を作りましょう．パラメータ`p`を引数に，とりうる値と確率関数を返します．
```python:ベルヌーイ分布の確率変数を返す関数
def Bern(p):
    x_set = np.array([0, 1])

    def f(x):
        if x in x_set:
            return (p ** x) * ((1 - p) ** (1 - x))
        else:
            return 0
    
    return x_set, f
```
関数が入れ子になっていて少し分かりにくいですが，ベルヌーイ分布の確率関数である`f(x)`を定義し，とりうる値と共に返しています．確率関数の実装は先程の`prob_bernoulli`と同じものです．同様に`p=0.2`（成功確率`20`%）で，確率変数の性質を確認する関数を実行してみます．
```pytbon:確率変数の性質の確認
p = 0.2
X = Bern(p)

check_prob(X)
```
```:実行結果
期待値: 0.2, 分散: 0.16
```
期待値は成功時は`1`，失敗時`0`になるので，`0.2 * 1 + 0.8 * 0 = 0.2`，分散は`0.2 * 0.8 = 0.16`になります．続いて，確率変数を図示します．

```python:確率変数の図示
plot_prob(X, size=(4, 4))
```
![スクリーンショット 2022-09-12 8.01.50.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/44273/cfeb5de9-bb8a-7086-8934-780a396ac3c3.png)

ベルヌーイ分布とその期待値の両方が表示できています．

### SciPyによる実装
続いて`scipy.stats`を使用した実装です．`scipy.stats`には，`p`を引数にベルヌーイ分布に従う確率変数に相当するオブジェクトを返す`bernoulli`関数が既に用意されています．

```python:SciPyでベルヌーイ分布に従う確率変数を返す
rv = stats.bernoulli(p)
rv
```
```:実行結果
<scipy.stats._distn_infrastructure.rv_frozen at 0x7f913611e990>
```
この`rv_frozen`は便利なメソッドを持っています．確率関数（Probability mass function）は`pmf`で，引数にとりうる値を指定すれば確率が返ってきます．

```python:確率関数
rv.pmf(0), rv.pmf(1)
```
```:実行結果
(0.8, 0.2)
```

とりうる値をリストで渡すこともできます．
```python:リストで渡す
rv.pmf([0, 1])
```
```:実行結果
array([0.8, 0.2])
```

累積分布関数（Cumulative distribution function）は`cdf`．
```python:累積分布関数
rv.cdf([0, 1])
```
```:実行結果
array([0.8, 0.2])
```

期待値は`mean`分散は`var`．
```python:期待値と分散
print(f'期待値:{rv.mean():.3f} 分散:{rv.var():.3f}')
```
```:実行結果
期待値:0.200 分散:0.160
```

### 使いどころ

ベルヌーイ分布に従う確率変数は，以下のように様々なケースに当てはめることができます．
- コイン投げで表が出るか
- サイコロを投げて6が出るか
- 広告がクリックされるか
- 来訪者が会員登録を完了するか
- LINEを送って返信が来るか

試行によって得られたデータから簡単に平均（期待値）が計算できるので，我々はベルヌーイ分布を特に意識することなく活用しているのかも，，，加えて，これから紹介する二項分布や幾何分布は，ベルヌーイ試行を何回か繰り返したときに現れる確率分布です．シンプルである分応用の幅が広い，確率分布の基本の1つと言えそうです．

ベルヌーイ分布の詳細については以下

https://www.headboost.jp/what-is-bernoulli/

## 二項分布（Binomial distribution）
成功確率`p`のベルヌーイ試行を`n`回行ったときの成功回数が従う分布を，**二項分布**と言います．パラメータはベルヌーイ分布と同じ`p`と，試行回数`n`です．とりうる値である成功回数は，1回も成功しない`0`から，すべて成功する`n`までのいずれかです．ちなみに，試行回数`n`が`1`のときの二項分布は，ベルヌーイ分布そのものになります！

例えば，表が出る確率が`1/2`のコインを`3`回投げたとき，とりうる値（成功回数）は以下になります．
- 成功回数: `0`・・・ぜんぶ裏が出る
- 成功回数: `1`・・・`1`回表が出て`2`回裏が出る
- 成功回数: `2`・・・`2`回表が出て`1`回裏が出る
- 成功回数: `3`・・・ぜんぶ表が出る

これらとりうる値それぞれの確率が，二項分布に従います．ベルヌーイ分布のときと同様に，まず確率変数を返す関数`Bin(n, p)`を作ります．

### NumPyによる実装
```python:二項分布の確率変数を返す関数
from scipy.special import comb

def Bin(n, p):
    x_set = np.arange(n + 1)
    
    def f(x):
        if x in x_set:
            return comb(n, x) * (p ** x) * ((1 - p) ** (n - x))
        else:
            return 0
        return x_set, f
    
    return x_set, f
```

ベルヌーイ分布の確率関数に似ているので，異なる部分を順番に説明していきます．

ベルヌーイ分布では，とりうる値は`[0, 1]`でしたが，以下のように`[0, 1, 2, ... n]`になります．
```python:ベルヌーイ分布のとりうる値
x_set = np.array([0, 1])
```
```python:二項分布のとりうる値
x_set = np.arange(n + 1)
```

確率関数の部分では，`scipy.special.comb`という関数を掛け合わせている点と，ベルヌーイ分布では累乗が`0`と`1`しかあり得なかったのが，上記のとりうる値の範囲で累乗される数式になっています．ベルヌーイ分布で妙に面倒な計算をしていたのは，二項分布の数式を基にしたからだったんですね，，，
```python:ベルヌーイ分布の確率関数
(p ** x) * ((1 - p) ** (1 - x))
```
```python:二項分布の確率関数
comb(n, x) * (p ** x) * ((1 - p) ** (n - x))
```
後半部分は，中学校で習った**積の法則**です．先程のコインを`3`回投げたときの例で考えると，成功回数が`1`のときは表が`1`回，裏が`2`回出る確率の積なので，`1/2 * 1/2 * 1/2 = 1/8`になります．成功回数のとりうる値すべてについて計算すると，以下のようになります．
- 成功回数`x = 0`・・・`(1/2 ** 0) * ((1 - 1/2) ** 3) = 1/8`
- 成功回数`x = 1`・・・`(1/2 ** 1) * ((1 - 1/2) ** 2) = 1/8`
- 成功回数`x = 2`・・・`(1/2 ** 2) * ((1 - 1/2) ** 1) = 1/8`
- 成功回数`x = 3`・・・`(1/2 ** 3) * ((1 - 1/2) ** 0) = 1/8`

コインは表と裏が出る確率が同じなので，すべて`1/8`になります．

前半部分も中学校で習った，**組み合わせ**と**和の法則**です．`3`回投げて成功回数が`1`のときに考えられる組み合わせは，「表，裏，裏」，「裏，表，裏」，「裏，裏，表」の`3`つなので，確率は`1/8 + 1/8 + 1/8 = 3/8`になります．`scipy.special.comb`はこの組み合わせを計算してくれる関数です．高校で習うコンビネーション（Combination），${}_nC_r$という奴ですね．とりうる値の確率をすべて計算すると以下のようになります．
- 成功回数`x = 0`・・・`comb(3, 0) * 1/8 = 1 * 1/8 = 1/8`
- 成功回数`x = 1`・・・`comb(3, 1) * 1/8 = 3 * 1/8 = 3/8`
- 成功回数`x = 2`・・・`comb(3, 2) * 1/8 = 3 * 1/8 = 3/8`
- 成功回数`x = 3`・・・`comb(3, 3) * 1/8 = 1 * 1/8 = 1/8`

積の法則，和の法則，組み合わせについては以下．

https://univ-juken.com/sekinohosoku-wanohosoku#i-2

https://blog.apar.jp/data-analysis/3927/

懐かしいですね．当時はこんなもの何の役に立つんだ..と文句ばかり言ってましたが，最近とても役立っています，，，当時の数学の先生ごめんなさい..

`Bin(n, p)`を使用してコインの例に従う確率変数を作ってみます．
```python:コインを3回投げたとき表が出る回数の確率分布
n = 3
p = 0.5
X = Bin(n, p)

x_set, bin = X
pd.DataFrame(
    [bin(x) for x in x_set],
    index=pd.Index(x_set, name='表が出る回数'),
    columns=['確率']).T
```
|表が出た回数|0|1|2|3|
|---|--:|--:|--:|--:|
|確率|0\.125|0\.375|0\.375|0\.125|

`0.125`は`1/8`，`0.375`は`3/8`なので，想定通りの確率分布を作ることができました．同様にサイコロを`10`回投げて`6`が出る回数の確率変数を作成してみます．

```python:サイコロを10回投げたとき6が出る回数の確率分布
n = 10
p = 1/6
X = Bin(n, p)

x_set, bin = X
prob = pd.DataFrame(
    [bin(x) for x in x_set],
    index=pd.Index(x_set, name='6が出る回数'),
    columns=['確率'])
prob.T
```
|6が出る回数|0|1|2|3|4|5|6|7|8|9|10|
|---|--:|--:|--:|--:|--:|--:|--:|--:|--:|--:|--:|
|確率|0\.162|0\.323|0\.291|0\.155|0\.054|0\.013|0\.002|0\.000|0\.000|0.000|0.000|

期待値と分散は以下のようになります．

```python:期待値と分散
check_prob(X)
```
```:実行結果
期待値: 1.667, 分散: 1.389
```
期待値は`n * p = 10 * 1 / 6 = 1.667`で，平均の成功回数を表します．これらの統計値はサイコロの出目の期待値や分散ではなく，`6`が出る**回数**の統計値である点に注意してください！確率分布と期待値を図示してみます．

```python:確率分布と期待値の図示
plot_prob(X)
```

![スクリーンショット 2022-10-02 12.00.33.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/44273/d4cb4083-f27e-1f15-e657-07e529b6d9d3.png)

サイコロを`10`回投げたとき`6`が出る回数の確率分布と期待値が確認できます．

例えば，`3`回成功する確率は約`15%`ということが分かります．また，`3`回以上成功する確率は，以下のように`3`回以上の確率を足し合わせれば算出できます．

```python:3回以上成功する確率
prob[prob.index >= 3].sum()
```
```:実行結果
確率   0.225
dtype: float64
```
`3`回以上成功する確率は，`22.5%`です！

### SciPyによる実装
SciPyでは`scipy.stats．binom`関数を使用するだけで二項分布の確率変数を得ることができます！`n`を`10`としたまま，`p`を`0.0〜1.0`で変化させて二項分布の形状を確認してみましょう．

```python:
n = 10
x_set = np.arange(n + 1)

plt.figure(figsize=(10, 6))
for i in range(1, 10):
    p = i / 10
    
    rv = stats.binom(n, p)
    prob = rv.pmf(x_set)

    plt.plot(x_set, prob, label=f'p:{p}')
plt.grid(True)
plt.xlabel('成功回数')
plt.ylabel('確率')
plt.legend()
plt.show()
```
![スクリーンショット 2022-10-02 17.31.29.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/44273/9c995529-387e-88f0-5c03-2bf04f4fe90e.png)

綺麗に左右対称になっていますね．`p`は成功確率なので，小さくなればなるほど`10`回試行したときの成功回数が少ない確率が高くなります．逆に`p`が大きくなればなるほど成功回数が多い確率が高くなります．

### 使いどころ
ベルヌーイ試行を，複数回行ったときの成功回数の確率分布が二項分布に従う，ということは，ベルヌーイ分布の使いどころで挙げたケースにおいて，成功回数の確率や期待値を求めたいときに役立ちそうです．例えば以下のようなケースを考えることができます．

内定率`15%`の学生が`10`社受けたとき，内定が`2`社以上得られる確率は？

まず，確率分布や期待値を求めます．
```python:内定社数の確率分布
X = Bin(0.15, 10)

check_prob(X)
plot_prob(X)
```
```:実行結果
期待値: 1.5, 分散: 1.275
```
![スクリーンショット 2022-10-10 12.10.01.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/44273/07e3fbf1-6039-373b-3a93-ed6c5e4d5676.png)
平均で`1.5`社の内定が期待できるということです．次に内定が`2`社以上得られる確率を出してみます．
```python:内定が2社以上得られる確率
x_set, bin = X
prob = [bin(x) for x in x_set[x_set >= 2]]
np.sum(prob)
```
```:実行時間
0.45570017623447245
```
内定が`1`社以下の値を`x_set`から除外した上で，`bin`関数に入力しています．約`46%`なりました！複数の内定から自分が行きたい会社を選びたいのであれば，もっと沢山の採用試験を受けたほうが良いのでは？と考えることができます．二項分布の詳細については以下．

https://www.headboost.jp/what-is-binomial-distribution/#index_id10

## 幾何分布（Geometric distribution）
ベルヌーイ試行を繰り返し，初めて成功したときの試行回数が従う確率分布が，**幾何分布**です．二項分布は一定回数試行したときの成功回数が従う確率分布ですが，幾何分布は一度成功したら終了します．逆に，成功しない限り制限無くベルヌーイ試行を繰り返します．従って，とりうる値は`1`以上の整数，パラメータは成功確率`p`だけです．

### NumPyによる実装
まず，NumPyで確率変数を実装します．無限のとりうる値を定義できないため，今回は上限の引数を設けて，デフォルト値を`100`としました．

```python:幾何分布の確率変数
def Ge(p, x_max=100):
    x_set = np.arange(1, x_max + 1)
    
    def f(x):
        if x in x_set:
            return p * ((1 - p) ** (x - 1))
        else:
            return 0
    return x_set, f
```
幾何分布についても，確率の積の法則に従った簡単な数式です．以下の部分です．
```python:
p * ((1 - p) ** (x - 1))
```
1番左の`p`は成功する確率です．`1`回成功したら終わりますが，成功するまで試行を続けるので，必ず`1`回は必要です．残りの`(1 - p) ** (x - 1)`は失敗する確率です．`1`回の成功以外はすべて失敗する確率なので，失敗する確率`1 - p`の累乗になります．これら2つの積を`x`回で成功する確率として返します．

`p`が`0.1`，`x`が`1〜3`のときの確率の計算は以下
- `1`回目で成功・・・`0.1 * (0.9 ** 0) = 0.1 * 1.0 = 0.1`
- `2`回目で成功・・・`0.1 * (0.9 ** 1) = 0.1 * 0.9 = 0.09`
- `3`回目で成功・・・`0.1 * (0.9 ** 2) = 0.1 * 0.81 = 0.081`

`2`回目で成功するということは，`1`回成功する確率と`1`回失敗する確率の積，`3`回目で成功する場合は，`1`回成功する確率と`2`回失敗する確率の積になります．回数を重ねるごとに徐々に確率が下がっていくことが分かります．

では，確率変数を作成してみます．
```python:エラーになる引数..
X = Ge(0.1, x_max=100)
check_prob(X)
```
```:実行結果
---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
<ipython-input-49-7d57ee186295> in <module>
      1 X = Ge(0.1, x_max=100)
----> 2 check_prob(X)

<ipython-input-5-fc7bec9aca03> in check_prob(X)
      6 
      7     prob_sum = np.round(np.sum(prob), 6)
----> 8     assert prob_sum == 1, f'確率の和が{prob_sum}になりました'
      9 
     10     print(f'期待値: {E(X):.4}, 分散: {V(X):.4}')

AssertionError: 確率の和が0.999973になりました
```
`check_prob`でエラーになってしまいました..とりうる値を`1〜100`に制限したため，確率の総和が`1.0`にならなかった模様..`x_max`を`1000`などと十分に大きくすれば問題無く動作しますが，今回は`p`を`0.2`とすることで対処しました．`p`を大きくすると，回数を重ねたときの確率がより小さくなり，合計が`1.0`に近づくため，エラーが起こらなくなりました．
```python:幾何分布の性質
X = Ge(0.2, x_max=100)
check_prob(X)
plot_prob(X, size=(24, 4))
```
```:実行結果
期待値: 5.0, 分散: 20.0
```
成功確率`20%`の試行を成功するまで続けたとき，平均すると`5`回で成功します．つまり，期待値は`1 / p`です．
![スクリーンショット 2022-10-11 10.50.24.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/44273/85a0443a-fa16-c62b-e16a-0aa32fa5ec74.png)
とりうる値を`1〜100`としたため，めっちゃ横に長いグラフになってしまった，，，試行が多いほど確率が小さくなっていくことがよく分かりますね．

### SciPyによる実装
`scipy.stats`による実装では，二項分布のときと同様に`p`を変化させて図示してみます．

```python:SciPyによる幾何分布
x_set = np.arange(21)

plt.figure(figsize=(10, 6))
for i in range(1, 10):
    p = i / 10
    
    rv = stats.geom(p)
    prob = rv.pmf(x_set)

    plt.plot(x_set, prob, label=f'p:{p}')
plt.grid(True)
plt.xlabel('試行回数')
plt.ylabel('確率')
plt.legend()
plt.show()
```
![スクリーンショット 2022-10-11 15.28.00.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/44273/04e3dfd5-3c3e-5036-1b2b-e265f5bb1983.png)
こちらは試行回数が少なくてもエラーにはなりませんでした，，，`p`が小さければなだらかに，`p`が大きければ急激に変化しています．成功確率が低ければ沢山試行しないと成功しない，成功確率が高ければすぐ成功して終わるという感覚と一致しますね！

### 使いどころ
幾何分布もベルヌーイ試行を繰り返しますが，一度でも成功したらそこで試行を終えます．例えば以下のようなケースを考えることができます．

1. サイコロを`1`が出るまで振ったとき，平均何回で終わるか
1. 当たるまで引ける当選確率`2`%の宝くじを`10`回以内に当てる確率は？
1. 継続率が`87`%の定期購入サービスの平均継続月数

1.は容易にイメージできるかと思います．2.は，`10`回までの確率分布の和を求めれば良さそうです．

3.の例を少し掘り下げてみます．継続率が`87`%ということは，解約率は`13`%です．ユーザーが解約することを一度決心し実行に移したら，そのユーザーに対する定期購入の売上はそこで途絶えてしまいます..つまり，この解約率を幾何分布のパラメータである`p`とすれば，成功するまでの試行回数が継続月数，期待値（平均）が平均継続月数になるはず！

では，平均継続月数の算出と，解約の確率分布の図示を行います．
```python:定期購入サービスの解約確率
x_set = np.arange(1, 32)

p = 0.13
rv = stats.geom(p)
pmf = rv.pmf(x_set)
e = rv.mean()
print(f'平均継続月数:{e:.03f}');

plt.figure(figsize=(10, 6))
plt.bar(x_set, pmf, label='解約確率')
plt.vlines(e, 0, 0.15, label='平均', color='red')
plt.grid(True)
plt.xticks(x_set)
plt.xlabel('継続月数')
plt.ylabel('確率')
plt.legend()
plt.show()
```
```:実行結果
平均継続月数:7.692
```
![スクリーンショット 2022-10-11 16.42.44.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/44273/bc628291-1cac-ae50-29c2-c8e405eb0260.png)
解約率`13`%のとき平均継続月数は約`8`ヶ月と求められ，解約の確率分布は上図のようになりました．赤い線として示した平均継続月数は，あくまで期待値（平均）として算出したものであって，解約確率は図のように分布していると考えることができます．実際の解約率`p`は固定ではなく，ユーザーや継続月数によって変化するケースが多いので，より正確な分布を得るには，もっと深く掘り下げて考える必要がありそうです，，，

この平均継続月数を使うと，例えば月額`1000`円のとき，LTV（Life time value:顧客生涯価値）を，`1000 * 7.692 = 7692`円と求めることができます．継続率のお話は以下のパチンコの例がわかり易かった，，，

https://note.com/mathpachi/n/nabd2f5bc0d27

幾何分布の詳細については以下

https://mathlandscape.com/geometric-distrib/#:~:text=%E5%B9%BE%E4%BD%95%E5%88%86%E5%B8%83%E3%81%A8%E3%81%AF%EF%BC%8C%E7%A2%BA%E7%8E%87,%E9%9B%A2%E6%95%A3%E5%9E%8B%E7%A2%BA%E7%8E%87%E5%A4%89%E6%95%B0%E3%81%A7%E3%81%99%E3%80%82

## ポアソン分布（Poisson distribution）
事象が単位時間あたりにランダムに発生するとき，発生件数が従う確率分布が**ポアソン分布**です．とりうる値は発生回数である`0`以上の整数です．引数は単位時間辺りの発生件数 λ（ラムダ）で正の実数，以降は`lam`とします．単位時間は`10`分でも`10`ヶ月でもなんでもOK．

二項分布のパラメータである試行回数`n`と成功確率`p`の積が平均成功回数ですが，これがポアソン分布での`lam`にあたります．言葉で表現しにくいのですが，一定期間常にベルヌーイ試行している状態で，何件か成功が発生するイメージ．つまり`n`が十分大きく，`p`が十分小さい状態です．このとき，`n * p`は一定期間辺りの発生件数として1つの数値で表すことができ，ポアソン分布ではこれを引数としているわけです！

### NumPyによる実装
では実装を行います．

```python:ポアソン分布の確率変数
from scipy.special import factorial

def Poi(lam, x_max):
    x_set = np.arange(x_max)
    
    def f(x):
        if x in x_set:
            return np.power(lam, x) / factorial(x) * np.exp(-lam)
        else:
            return 0
    return x_set, f 
```
難しい関数がいくつか出てきましたね..`np.power`は累乗を求める関数です．`scipy.special.factorial`は階乗を求める関数です．そして`np.exp`はネイピア数の累乗です．ここでは式の細かい説明は割愛しますが，二項分布からポアソン分布を導出することができます．詳細は以下．

https://note.com/konpyu/n/ncede5754aba8

確率変数を作成して性質を確認します．
```python:ポアソン分布の性質
lam = 4
X = Poi(4, x_max = 30)

check_prob(X)
```
```:実行結果
期待値: 4.0, 分散: 4.0
```
二項分布と同様に期待値は`n * p`なので，そのまま`lam`の値になります．確率分布を図示してみましょう．

```python:ポアソン分布の図示
plot_prob(X)
```
![スクリーンショット 2022-10-13 14.58.33.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/44273/20ee1319-482f-77f0-830d-87a0109ec14b.png)

### SciPyによる実装
`scipy.stats.poisson`でポアソン分布の確率変数を作ることができます．

```python:lamを変化させたときのポアソン分布
x_set = np.arange(31)

plt.figure(figsize=(10, 6))
for lam in range(1, 20, 2):
    rv = stats.poisson(lam)
    prob = rv.pmf(x_set)

    plt.plot(x_set, prob, label=f'lam:{lam}')
plt.grid(True)
plt.xlabel('発生件数')
plt.ylabel('発生確率')
plt.legend()
plt.show()
```
![スクリーンショット 2022-10-13 15.06.36.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/44273/834e6a10-4769-766f-09f4-f84923489408.png)
`lam`を`1〜19`で変化させて図示しました．単位時間辺りの発生件数`lam`を大きくするほど分布が広がっています．また，確率の最大値が`lam`に一致します．例えば以下のようなケースが考えられます．

### 使いどころ
ポアソン分布は，**試行のタイミングが明確ではない**が，一定期間にランダムに発生する事象について，その発生回数が従う確率分布です．例えば以下のようなケースを考えることができます．ただし，一定期間辺りの発生件数には平均を用いているので，時間帯やトレンドを考慮したい場合は，発生件数の変化を考慮に入れる必要があります．

1. 年間`0.61`人が馬に蹴られて死亡する軍隊において，年間何人が馬に蹴られて死ぬか
1. `15`分の間に平均`8`件電話がかかってくるコールセンターにて，`15`分間に`17`件以上かかってくる確率


1.は歴史的に有名な事例らしいです．図示してみました．

```python:年間で落馬で死亡する人数の確率分布
x_set = np.arange(11)
lam = 0.61

rv = stats.poisson(lam)
prob = rv.pmf(x_set)

plt.figure(figsize=(10, 6))
plt.bar(x_set, prob)
plt.grid(True)
plt.xlabel('死亡件数')
plt.ylabel('発生確率')
plt.show()
```
![スクリーンショット 2022-10-13 15.41.56.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/44273/c5cf73c1-9684-8baf-96af-1d5007d88638.png)
平均発生件数が`1.0`未満のケースでも，問題無く確率分布を出すことができていますね．

2.は現代でも役立ちそうなケース．

```python:コールセンターのキャパシティ
x_set = np.arange(40)
lam = 8

rv = stats.poisson(lam)
prob = rv.pmf(x_set)

plt.figure(figsize=(10, 6))
plt.bar(x_set, prob)
plt.grid(True)
plt.xlabel('コール件数')
plt.ylabel('発生確率')
plt.show()
```
![スクリーンショット 2022-10-13 15.48.58.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/44273/45397306-5ba5-3fef-d838-27177b41f8a3.png)

`17`件以上の確率を合計してみます．
```python:コールセンターのキャパシティの検証
np.sum(prob[17:])
```
```:実行結果
0.0037180212812834303
```
例えば１人辺り`15`分で`8`件電話応対できる人が`2`名いたとき，電話応対のキャパシティは`16`件です．この場合，15分辺り約`0.37`%の確率でキャパシティをオーバーする可能性がある，逆に`99.63`%はキャパシティ内に収まる，ということが分かります．

電話の向こうにいる大勢の人が，何をどのタイミングで何回試行しているのか，試行してどんな問題が発生したのか，電話がかかってこない限り，コールセンター側は知ることができません．しかし，一定時間に何件電話がかかってきたか，という情報は容易に得ることができます．ポアソン分布を使うことで，その情報から上記のような分析が可能になるのです！

他にも機器の故障や，多チャネルで販売されている商品の購入など，試行回数の把握が難しいケースで活用されているようです．ポアソン分布の詳細は以下

https://atarimae.biz/archives/7372

# 連続型確率変数
ここまで離散型の確率変数を扱ってきましたが，ここからは，とりうる値が連続的な確率変数について説明します．離散型の確率変数は，とりうる値が飛び飛びの点なので，対応する確率分布も飛び飛びの点になります．一方，連続型の確率変数では，とりうる値が連続した線なので，対応する確率も連続した線になります．

## 面積と積分
例えばサイコロであれば，`2`か`3`が出る確率は，`1/6 + 1/6`で`1/3`で，飛び飛びの確率を足せば求めることができます．一方で，連続型の確率変数では，とりうる値が連続しているので，`2`と`3`の間に`2.5`や`2.999`といったとりうる値が無限に存在するため，とりうる値に対応する確率を足して求めることができません，，，

そこで，`2〜3`というように，とりうる値の範囲を指定して，その部分の確率分布の面積を求めることで，とりうる値が指定した範囲に収まる確率を算出することができます！面積を求める手法として，高校のときに習った**積分（Integral**）を用いています．積分に関する数学的な説明は省略します．

積分というワードを敬遠する方がいるかもしれませんが，連続型確率分布において，とりうる値の範囲の面積を計算することで確率を求めることができる，ということが理解できれば，問題無く読み進めることができるかと思います．

積分に関する簡単な説明は以下

https://www.headboost.jp/what-is-integral/

## 1次元
連続型確率変数の説明は，確率編1で説明した[離散型確率変数](https://qiita.com/o93/items/1b3f379da3e62c6252d8#%E9%9B%A2%E6%95%A3%E5%9E%8B%E7%A2%BA%E7%8E%87%E5%A4%89%E6%95%B0discrete-random-variable)と同じ構成になっています．積分を用いたり，使用する計算式は離散型と違いますが，確率分布の考え方は同じですので，離散型と見比べながら確認すると理解しやすいかもしれません．

1次元の連続型確率変数では具体例として**連続型いかさまルーレット**を考えることにします．値を指し示す針があり，ルーレットをくるくる回します．円は1周の長さを`1.0`とし，ルーレットが止まったときの始点から針までの弧の長さを出目（実現値）とします．つまり，ルーレットのとりうる値は`0.0〜1.0`の間の実数になります．そして，離散型確率変数の際に扱った「いかさまサイコロ」と同様に，大きい値ほど出やすい，といういかさま仕様です，，，

以下のようなイメージ．
![スクリーンショット 2022-10-23 8.26.33.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/44273/a010541d-69ed-9251-80ed-97339e11cb32.png)

さて，ルーレットを回して↑のようになる確率を考えてみましょう．`0.58..`辺りでしょうか，，，通常の離散型ルーレットのように出目がはっきり区切られていれば，値を定めることができますが，連続型ルーレットには区切りが存在しないため，値が無限に存在します．無限に存在する値の確率が`0`より大きいと，確率の合計が`1.0`になるという，確率変数の条件を満たすことができません..仮にとりうる1つの値を指定すると，確率は`0`になってしまいます．

そこで連続型確率変数では，とりうる値を区間として指定し，その区間に入る確率を考えます．例えば，ルーレットが`0.5〜0.7`を指す確率は，`0.3`といった感じです．

### 確率密度関数／密度関数（Probability denssity function／PDF）
離散型確率変数は，とりうる値と確率質量関数（確率関数）で定義できましたが，連続型確率変数は，とりうる値の区間と**確率密度関数（密度関数**）で定義します．ここで重要なのは，確率密度関数と確率質量関数の違いです．確率密度関数はとりうる値の確率を返すわけではなく，とりうる値の**確率密度**を返します！

これを理解するには，関数の名前にもあるように，**質量**と**密度**の違いを考える必要があります．中学生の頃に習った公式を思い出してみましょう．

`体積 * 密度 = 質量`

同じ体積の鉄と木材では鉄のほうが重い．何故なら，鉄の密度が高いからです．質量は体積と密度を掛け合わせて求めることができます．質量と体積の関係については以下．

http://kentiku-kouzou.jp/suugaku-mitudotaiseki.html

確率もこれと同じように考えると，以下のようになります．

`とりうる値の区間 * 確率密度 = 確率`

とりうる値の区間と確率密度をかけ合わせることで，とりうる値がその区間内に収まる確率を求めることができます．離散型はとりうる値が飛び飛びなので，区間がありません．なので，値から直接確率を算出することができます．一方で連続型は値が連続しているため，値から直接確率を算出しようとすると，ルーレットの例で挙げたようにすべて`0`になってしまいます．そこで，とりうる値の相対的な出やすさである確率密度を返す関数を定義し，確率密度ととりうる値の区間の面積を確率として表現しているのです！

ここからは，いかさまルーレットを例にして密度関数を実装していきます．まず，とりうる値の区間を定義します．
```python:いかさまルーレットのとりうる値の区間
x_range = np.array([0, 1])
```
最低値と最大値を指定しています．

次に密度関数を定義します．
```python:いかさまルーレットの密度関数
def f(x):
    x_min, x_max = x_range
    if x_min <= x <= x_max:
        return 2 * x
    else:
        return 0
```
`x`が`x_range`の範囲内であるとき，`2 * x`を返し，それ以外は`0`を返します．いかさまサイコロのときと同様に，`x`が増えるほど大きい確率を返す仕組みになっていますが，`2`をかけているのは，確率の性質である，全確率の総和を`1.0`にするための工夫です．

では，確率変数を定義して密度関数を図示してみましょう！
```python:密度関数と出目が0.5〜0.7になる確率
X = [x_range, f]

x_min, x_max = x_range
x_s = np.linspace(x_min, x_max, 100)

x_fill = np.linspace(0.5, 0.7, 100)

plt.figure(figsize=(10, 6))
plt.plot(x_s, [f(x) for x in x_s], label='確率密度関数')
plt.fill_between(x_fill, [f(x) for x in x_fill], label='確率', alpha=0.3)

plt.grid(True)
plt.xlabel('とりうる値')
plt.ylabel('密度')
plt.legend()
plt.show()
```
![スクリーンショット 2022-10-23 10.53.59.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/44273/14a384a5-4910-c55c-abb5-9c75a7c0d964.png)
`x_s`や`x_fill`は，グラフを描画するためだけに使用する配列です．いかさまルーレットの出目が`0.5〜0.7`の範囲内に収まる確率の領域を塗りつぶしています．離散型では縦軸が確率になっていましたが，連続型では，縦軸は**確率密度**を表しており，塗りつぶされた領域の面積が確率になっている点に注意してください！

とりうる値と確率密度が得られたので，塗りつぶした領域の面積（確率）を求めていきましょう．積分には`scipy.integrate.quad`を使用します．

```python:いかさまルーレットの確率分布の面積
integrate.quad(f, 0.5, 0.7)
```
```:実行結果
(0.240, 0.000)
```
積分値と推定誤差を返してくれます．確率は`0.24`と出ました．台形の面積を求める公式を使って，積分によって得られた確率が正しいか，確認してみます．区間は`0.5〜0.7`なので，台形の高さは`0.2`，上底は`1.0`，下底は`1.4`なので，

`(1.0 + 1.4) * 0.2 * 0.5 = 2.4 * 0.1 = 0.24`
おお〜積分の結果と台形の面積が一致しましたね！

確率の性質は離散型と同じで，すべて`0.0`以上，総和は`1.0`になります．
```python:確率の性質
f_min = optimize.minimize_scalar(f).fun
f_sum = integrate.quad(f, -np.inf, np.inf)[0]

print(f'最小値: {f_min}')
print(f'総和: {f_sum:.3f}')
```
```:実行結果
最小値: 0
総和: 1.000
```
`scipy.optimize.minimize_scalar`は関数の戻り値の最小値を求める関数で，確率密度関数の最小値が`0`になることを求めています．総和は先程使用した`scipy.integrate.quad`にて，区間を無限大にした積分値を求めることができます．

総和についても，今回のケースだと三角形の面積で求めることができます．とりうる値の区間は`1.0`，確率密度の高さは`2.0`なので，

`1.0 * 2.0 * 0.5 = 1.0`
となり，総和と一致します．例えば密度関数を単に`x`とした場合，三角形の面積は`0.5`になってしまい，確率の性質に一致しません，，，確率の性質に一致させるために，密度関数内で`2 * x`としていたわけですね！

### 累積分布関数／分布関数（Cumulative distribution function／CDF）
**累積分布関数**は，確率変数Xが引数x以下になる確率を返します．こちらは離散型確率変数のときと同じです．

```python:累積分布関数
def F(x):
    return integrate.quad(f, -np.inf, x)[0]

F(0.7)
```
```:実行結果
0.490
```
先程の`scipy.integrate.quad`を使用して実装しています．例えば，いかさまルーレットが`0.7`以下になる確率は`49`%と求めることができます．分布関数を組み合わせることで，区間の指定も可能です．

```python:累積分布関数による区間指定
F(0.7) - F(0.5)
```
```:実行結果
0.240
```

分布関数を図示してみます．
```python:いかさまルーレットの累積分布関数の図示
x_min, x_max = x_range
x_s = np.linspace(x_min, x_max, 100)

plt.figure(figsize=(10, 6))
plt.plot(x_s, [F(x) for x in x_s], label='累積分布関数')

plt.grid(True)
plt.xlabel('とりうる値')
plt.ylabel('確率')
plt.legend()
plt.show()
```
![スクリーンショット 2022-10-23 12.17.21.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/44273/cfaf49e1-32ae-41c0-87c7-58bbefee54e2.png)
確率を累積しているため，`0`から始まって，とりうる値の最大値`1.0`で`1.0`になっています．確率密度関数と違って，こちらは縦軸が密度ではなく**確率**になっている点に注意してください．

### 期待値（Expected value）
連続型確率変数の**期待値**（平均）です．前記事の[離散型確率変数の期待値](https://qiita.com/o93/items/1b3f379da3e62c6252d8#%E6%9C%9F%E5%BE%85%E5%80%A4expected-value)と同じで，とりうる値と確率の積の総和になります．期待値を求める関数を実装します．

```python:確率の期待値を求める
def integrand(x):
    return x * f(x)

integrate.quad(integrand, -np.inf, np.inf)[0]
```
```:実行結果
0.667
```
いかさまルーレットの密度関数にとりうる値を掛けて，その全面積を求めることで総和を算出できます．これが期待値です．離散型と同じく，期待値を求める関数`E(X)`を作ります．
```python:期待値を求める関数
def E(X, g=lambda x: x):
    x_range, f = X

    def integrand(x):
        return g(x) * f(x)
    
    return integrate.quad(integrand, -np.inf, np.inf)[0]

E(X)
```
```:実行結果
0.667
```
[離散型確率変数の変換](https://qiita.com/o93/items/1b3f379da3e62c6252d8#%E7%A2%BA%E7%8E%87%E5%A4%89%E6%95%B0%E3%81%AE%E5%A4%89%E6%8F%9B)のように，確率変数の変換にも対応しています．例えば，いかさまルーレットの出目に対して以下のルールに基づいて報酬を払うことにします．

- 出目の`5`倍の報酬が貰える
- 参加費用として`3`取られる

```python:確率変数の変換
E(X, g=lambda x: 5 * x - 3)
```
```:実行結果
0.333
```
引数`g`に報酬関数を指定することで，期待値を算出することができます！こちらも離散型のときと同様に，線形性が成り立ちます．

```python:期待値の線形性
5 * E(X) - 3
```
```:実行結果
0.333
```

### 分散（Variance）
連続型確率変数の**分散**についても，式は[離散型の分散](https://qiita.com/o93/items/1b3f379da3e62c6252d8#%E5%88%86%E6%95%A3variance)と同じです．同様に`V(X)`を実装します．

```python:分散を求める関数
def V(X, g=lambda x: x):
    x_range, f = X
    mean = E(X, g)

    def integrand(x):
        return (g(x) - mean) ** 2 * f(x)
    
    return integrate.quad(integrand, -np.inf, np.inf)[0]

V(X)
```
```:実行結果
0.056
```

期待値のときと同じように`g`を指定してみます．
```python:確率変数の変換
V(X, g=lambda x: 5 * x - 3)
```
```:実行結果
1.389
```
また，離散型でも示した分散の公式は，連続型でも成り立ちます．
```python:分散の公式
5 ** 2 * V(X)
```
```:実行結果
1.389
```
同じ値になりましたね！

## 2次元
次は2次元の連続型確率変数について扱っていきます．

### 同時確率密度関数（Joint probability density function）
2つのとりうる値の区間の組み合わせの密度を返す関数が，**同時確率密度関数**です．離散型の2次元では2つのいかさまサイコロを扱いました．連続型でも２つのいかさまルーレットA，Bを使って考えます．ルーレットAの確率変数を`Y`，ルーレットAとルーレットBの値を足したものを確率変数`X`とします．[離散型のとき](https://qiita.com/o93/items/1b3f379da3e62c6252d8#%E5%90%8C%E6%99%82%E7%A2%BA%E7%8E%87%E5%88%86%E5%B8%83%E5%90%8C%E6%99%82%E5%88%86%E5%B8%83joint-probability-distribution)と同じように，あとあと共分散や相関係数を求める都合で，２つの変数に関連をもたせています．まとめると以下のようになります．
- `X`・・・「ルーレットAの値 + ルーレットBの値」の確率変数
- `Y`・・・「ルーレットBの値」の確率変数

いかさまルーレットのとりうる値の範囲は`0.0〜1.0`なので，とりうる値`x`，`y`のそれぞれの範囲は以下．
- `x`のとりうる範囲・・・`0.0〜2.0`
- `y`のとりうる範囲・・・`0.0〜1.0`

ここで，同時確率密度関数の実装をどのようにすべきか考えます．書籍では何故か詳しく触れられておらず，いきなり実装の式と条件が出てきて混乱してしまいましたが，，，

最初にルーレットの値を確率変数のとりうる値`x`，`y`から逆算します．実装の都合上，いかさまルーレットA，Bそのものではなく，`x`と`y`から確率密度を導き出したいためです．`x = ルーレットAの値 + ルーレットBの値`なので，式変形して`ルーレットAの値 = x - ルーレットBの値`．更に，`ルーレットBの値 = y`なので，以下のようにルーレットA，Bの値を`x`と`y`で表すことができます．

- ルーレットAの値・・・`x - y`
- ルーレットBの値・・・`y`

いかさまルーレットの値の範囲は`0.0〜1.0`なので，
- ルーレットAのとりうる範囲・・・`0 <= x - y <= 1`
- ルーレットBのとりうる範囲・・・`0 <= y <= 1`

になります．2つのいかさまルーレットの確率密度は1次元のときの例と同じで，いかさまルーレットの値の2倍です．
- ルーレットAの確率変数・・・`2 * (x - y)`
- ルーレットBの確率変数・・・`2 * y`

確率の積の法則から，
`2 * (x - y) * 2 * y = 4 * y * (x - y)`

従って，確率変数`(X, Y)`の同時確率密度関数は以下のようになります！
- `0 <= y <= 1`かつ`0 <= x - y <= 1`のとき・・・`4 * y * (x - y)`
- それ以外のとき・・・`0`

上記を実装してみましょう．

```python:とりうる値の定義
x_range = [0, 2]
y_range = [0, 1]
```
```python:同時確率密度関数の定義
def f_xy(x, y):
    if 0 <= y <= 1 and 0 <= x - y <= 1:
        return 4 * y * (x - y)
    else:
        return 0
```
```python:2次元の確率変数
XY = [x_range, y_range, f_xy]
XY
```
```:実行結果
[[0, 2], [0, 1], <function __main__.f_xy(x, y)>]
```
同時確率密度関数を使って，確率密度の分布をヒートマップで描画します．

```python:確率密度の分布
xr = np.linspace(x_range[0], x_range[1], 200)
yr = np.linspace(y_range[0], y_range[1], 200)
pd_xy = np.array([[f_xy(x, y) for y in yr] for x in xr])

plt.figure(figsize=(10, 8))
sns.heatmap(
    pd_xy, xticklabels=np.linspace(0, 2, 3), yticklabels=np.linspace(0, 1, 3),
    cmap='Blues')
plt.xticks(np.linspace(0, 200, 3))
plt.yticks(np.linspace(0, 200, 3))
plt.xlabel('X')
plt.ylabel('Y')
plt.show()
```
![スクリーンショット 2022-10-29 11.44.06.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/44273/a5b95f79-c6ac-931b-2e11-ca0c6c267d4d.png)
色の濃さが確率の密度を表しています．確率ではなく確率の**密度**であることに注意してください．右側にあるバーから，確率が必ず`0`以上になることが分かります．ちなみに[離散型のヒートマップ](https://camo.qiitausercontent.com/f1f08797cbb4fa302d893611583e0fd319c4dd3b/68747470733a2f2f71696974612d696d6167652d73746f72652e73332e61702d6e6f727468656173742d312e616d617a6f6e6177732e636f6d2f302f34343237332f32356566346234352d376462322d316637322d636165632d3432616666656635303164652e706e67)に似ていますよね．これは，確率変数の性質が同じだからです．離散型も`x`はサイコロA + サイコロBの出目で，`y`はサイコロBの出目なので．

確率の総和は，`integrate.nquad`で求めます．`X`を横軸，`Y`を縦軸，確率密度を高さとすれば，確率は立体の体積に相当します．

```python:2次元の確率の総和
integrate.nquad(f_xy, [[-np.inf, np.inf], [-np.inf, np.inf]][0]
```
```:実行結果
1.000
```
以上で，確率の条件を満たしていることが確認できました．

### 周辺確率密度関数（Marginal probability density function）
[離散型の周辺分布](https://qiita.com/o93/items/1b3f379da3e62c6252d8#%E5%91%A8%E8%BE%BA%E7%A2%BA%E7%8E%87%E5%88%86%E5%B8%83%E5%91%A8%E8%BE%BA%E5%88%86%E5%B8%83marginal-probability-distribution)と同様に，複数の物事からなる確率密度関数から，1つの確率密度の変化だけを抽出した関数を，**周辺確率密度関数（周辺密度関数**）と言います．離散型のときと同様に片方の引数を固定して求めることができます．

```python:周辺密度関数の定義
def f_X(x):
    return integrate.quad(partial(f_xy, x=x), -np.info, np.inf)[0]

def f_Y(y):
    return integrate.quad(partial(f_xy, ｙ=y), -np.info, np.inf)[0]
```
`f_xy`の片方の引数を固定するため，標準ライブラリ`functools.partial`を使用しています．片方の引数を固定した確率密度関数を`integrate.quad`に指定して，その面積を求めています．それでは，`X`，`Y`それぞれの密度関数を図示してみましょう．

```python:XとYそれぞれの周辺密度関数の図示
X = [x_range, f_X]
Y = [y_range, f_Y]

xr = np.linspace(x_range[0], x_range[1], 100)
yr = np.linspace(y_range[0], y_range[1], 100)

prob_x = [f_X(x) for x in xr]
prob_y = [f_Y(y) for y in yr]

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

ax1.plot(xr, prob_x)
ax1.set_title('Xの密度関数')
ax1.set_xlabel('Xのとりうる値')
ax1.set_ylabel('密度')
ax1.grid()

ax2.plot(yr, prob_y)
ax2.set_title('Yの密度関数')
ax2.set_xlabel('Yのとりうる値')
ax2.set_ylabel('密度')
ax2.grid()

plt.show()
```
![スクリーンショット 2022-10-29 12.53.49.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/44273/af6ab7d6-0a30-dc00-7ce8-81d41099d33e.png)
こっちも[離散型のときの周辺分布](https://camo.qiitausercontent.com/6ed78e00fd267e477e31ebabf1e7f28b7e1ebbd6/68747470733a2f2f71696974612d696d6167652d73746f72652e73332e61702d6e6f727468656173742d312e616d617a6f6e6177732e636f6d2f302f34343237332f66396636396562642d623566362d613565312d343162302d3165643436323135303966322e706e67)に似ていますね！

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

```python:2次元の連続型確率変数の期待値
def E2(XY, g):
    x_range, y_range, f_xy = XY
    
    def integrand(x, y):
        return g(x, y) * f_xy(x, y)
    
    return integrate.nquad(integrand, [[-np.inf, np.inf], [-np.inf, np.inf]])[0]

mean_X = E2(XY, lambda x, y: x)
mean_Y = E2(XY, lambda x, y: y)
mean_X, mean_Y
```
```:実行結果
(1.333, 0.667)
```
[離散型の2次元の期待値](https://qiita.com/o93/items/1b3f379da3e62c6252d8#%E6%9C%9F%E5%BE%85%E5%80%A4expected-value-1)と似たようなコードですね．同様に，期待値の線形性も成り立ちます．

```python:期待値の線形性
a, b = 2, 3

E2(XY, lambda x, y: (a * x) + (b * y)), (a * mean_X) + (b * mean_Y)
```
```:実行結果
(4.667, 4.667)
```

### 分散（Variance）
分散も1次元のときとほぼ同じです．

```python:2次元の連続型確率変数の分散
def V2(XY, g):
    x_range, y_range, f_xy = XY
    mean = E2(XY, g)
    
    def integrand(x, y):
        return (g(x, y) - mean) ** 2 * f_xy(x, y)
    
    return integrate.nquad(integrand, [[-np.inf, np.inf], [-np.inf, np.inf]])[0]

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

### 共分散（Covariance）
共分散はXとYそれぞれの偏差の積の期待値によって求められ．2つの確率変数にどの程度相関があるかを示します．こちらも[離散型の2次元の共分散](https://qiita.com/o93/items/1b3f379da3e62c6252d8#%E5%85%B1%E5%88%86%E6%95%A3covariance)と同じ．

```python:2次元の連続型確率変数の共分散
def cov(XY):
    x_range, y_range, f_xy = XY
    mean_X = E2(XY, lambda x, y: x)
    mean_Y = E2(XY, lambda x, y: y)
    
    def integrand(x, y):
        return (x - mean_X) * (y - mean_Y) * f_xy(x, y)
    
    return integrate.nquad(integrand, [[-np.inf, np.inf], [-np.inf, np.inf]])[0]

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

### 相関係数（Correlation coefficient）
相関係数は，それぞれの分散と共分散から求めることができます．これも[離散型の2次元の相関係数](https://qiita.com/o93/items/1b3f379da3e62c6252d8#%E7%9B%B8%E9%96%A2%E4%BF%82%E6%95%B0correlation-coefficient)と同じです．

```python:2次元の連続型確率変数の相関係数
cov_XY / np.sqrt(var_X * var_Y)
```
```:python
0.707
```
離散型と同じ値になりましたね！！いかさまサイコロA，BといかさまルーレットA，Bの確率変数は本質は同じ確率分布に従っているため，離散，連続の違いはあれどAとBの相関係数も同様の値になるわけです．

# おわりに
ここまでを確率編2とします！今回は，代表的な離散型確率分布の紹介と使いどころ，連続型確率変数の基本について記載しました．

Androidエンジニアをやっていた頃は正規分布くらいしか知らなかったのですが，色んな確率分布があり，それぞれの仕組み・用途・関係性を改めて確認することができました．

また，中学・高校の頃に習った確率の和・積の法則や，組み合わせ（コンビネーション），密度と質量，面積と積分など，勉強嫌いだった当時は一生使うことは無いと思っていたものが，数十年後，業務に役立つ技術の理解を助けてくれるとは思いもせず，，，

しかし，記事を書いている自分自身，この考え方は何の役に立つんだ..とか，この数式になんの意味があるんだ..とか，？が沢山浮かんできます．その答えがすべて書籍に書かれているわけもなく，なるべく調べて補足を入れているものの，依然いくつかの疑問は残ったままです．引き続き勉強を続けて，「なるほど！」となったタイミングで記事を修正していければと思います．記事に不備や問題があればいつでもご連絡ください！

次回は確率編のおそらくラストで，代表的な連続型確率分布，多次元確率変数などについて記載する予定です．