見出し画像

StanとRでベイズ統計モデリングをPyMC Ver.5で写経～第６章「さまざまな確率分布」

ネイピア DS

第６章「統計モデリングの視点から確率分布の紹介」

書籍の著者　松浦健太郎 先生

この記事は、テキスト第６章「統計モデリングの視点から確率分布の紹介」の Python 写経 を取り扱います。
ベイズモデリングの大切な仲間【確率分布】の確率密度関数・確率質量関数を描きます。
Python の確率・統計ライブラリ scipy.stats を利用します。
三角図描画は先達のお知恵を拝借し、条件付き確率は著者先生のコードをお借りしました。

はじめに

StanとRでベイズ統計モデリングの紹介

この記事は書籍「StanとRでベイズ統計モデリング」（共立出版、「テキスト」と呼びます）のベイズモデルを用いて、PyMC Ver.5で「実験的」に写経する翻訳的ドキュメンタリーです。

テキストは、2016年10月に発売され、ベイズモデリングのモデル式とプログラミングに関する丁寧な解説とモデリングの改善ポイントを網羅するチュートリアル「実践解説書」です。もちろん素晴らしいです！
「アヒル本」の愛称で多くのベイジアンに愛されてきた書籍です！

テキストに従ってStanとRで実践する予定でしたが、RのStan環境を整えることができませんでした（泣）
そこでこのシリーズは、テキストのベイズモデルをPyMC Ver.5に書き換えて実践します。

引用表記

この記事は、出典に記載の書籍に掲載された文章及びコードを引用し、適宜、掲載文章とコードを改変して書いています。
【出典】
「StanとRでベイズ統計モデリング」初版第13刷、著者 松浦健太郎、共立出版

記事中のイラストは、「かわいいフリー素材集いらすとや」さんのイラストをお借りしています。
ありがとうございます！


PyMC環境の準備

Anacondaを用いる環境構築とGoogle ColaboratoryでPyMCを動かす方法について、次の記事にまとめています。
「PyMCを動かすまでの準備」章をご覧ください。


６章 さまざまな確率分布

インポート

### インポート

# 数値・確率計算
import pandas as pd
import numpy as np
import scipy.stats as stats

# 描画
import matplotlib.pyplot as plt
from matplotlib import gridspec
plt.rcParams['font.family'] = 'Meiryo'

# ワーニング表示の抑制
import warnings
warnings.simplefilter('ignore')

6.1 一様分布

図6.1を描画します。
確率分布には scipy.stats の uniform() を利用します。

### 図6.1の描画

## 設定
# パラメータa, b
a, b = -2, 1
# 線・マーカーの色
c = 'tab:blue'

## 描画用データの作成
# 一様分布のパラメータscaleの計算：a=loc, b=loc+scale, scale=b-a 
scale = b - a
# 確率変数yの値（確率密度が0以外の区間）
yvals = np.linspace(a, b, 1001)
# 上記の確率変数yの値に対する確率密度関数の算出
density = stats.uniform.pdf(x=yvals, loc=a, scale=scale)

## 描画処理
# 描画領域の設定
plt.figure(figsize=(5, 3))
ax = plt.subplot()
# -2から1の区間の描画 線・両端の点
ax.plot(yvals, density, color=c)
ax.plot(yvals[0], density[0], 'o', ms=8, color=c)
ax.plot(yvals[-1], density[-1], 'o', ms=8, color=c)
# -3から-2の区間の描画  線・片端の白抜き点
ax.plot([-3, a], [0, 0], color=c)
ax.plot(a, 0, 'o', ms=8, color=c, mfc='white')
# 1から3の区間の描画  線・片端の白抜き点
ax.plot([b, 3], [0, 0], color=c)
ax.plot(b, 0, 'o', ms=8, color=c, mfc='white')
# 修飾
ax.set(title=f'一様分布\n$a={a},\ b={b}$', xlabel='確率変数 $y$', ylabel='density',
       xlim=(-3, 3), ylim=(-0.05, 0.38))
plt.grid(lw=0.5);

【実行結果】

6.2 ベルヌーイ分布

図6.2を描画します。
確率分布には scipy.stats の bernoulli() を利用します。

### 図6.2の描画

## 設定
# パラメータ theta
theta = 0.2
# 線・マーカーの色
c = 'tab:blue'

## 描画用データの作成
# 確率変数yの値
yvals = [0, 1]
# 確率変数yの値に対する確率質量関数の算出
prob = stats.bernoulli.pmf(k=yvals, p=theta)

## 描画処理
# 描画領域の設定
plt.figure(figsize=(5, 3))
ax = plt.subplot()
# 棒グラフの描画
ax.bar(yvals, prob, color=c)
# 修飾
ax.set(title=f'ベルヌーイ分布\n$p={theta}$', xlabel='確率変数 $y$',
       ylabel='probability', xticks=(0, 1), ylim=(-0.05, 0.87))
plt.grid(lw=0.5);

【実行結果】

6.3 二項分布

図6.3を描画します。
確率分布には scipy.stats の binom() を利用します。

### 図6.3の描画

## 設定
# パラメータ
N, theta = 10, 0.2
# 線・マーカーの色
c = 'tab:blue'

## 描画用データの作成
# 確率変数yの値
yvals = range(0, 11)
# 確率変数yの値に対する確率質量関数の算出
prob = stats.binom.pmf(k=yvals, n=N, p=theta)

## 描画処理
# 描画領域の設定
plt.figure(figsize=(5, 3))
ax = plt.subplot()
# 棒グラフの描画
ax.bar(yvals, prob, color=c)
# 修飾
ax.set(title=f'二項分布\n$N={N},p={theta}$',
       xlabel='確率変数 $y$',ylabel='probability',
       xticks=yvals, ylim=(-0.01, 0.33), yticks=(0, 0.1, 0.2, 0.3))
plt.grid(lw=0.5);

【実行結果】

6.4 ベータ分布

図6.4を描画します。
確率分布には scipy.stats の beta() を利用します。

### 図6.4の描画

## 設定
# 色
c = ['tab:blue', 'tab:red', 'tab:green']
# 線種
ls = ['-', '--', '-.']

## 描画用データの作成
# 確率変数yの値
yvals = np.linspace(0, 1, 1001)

## 描画処理
# 描画領域の設定
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(9, 5))
# 左のグラフの描画
for i, (a, b) in enumerate(zip([1, 3, 3], [1, 2, 9])):
    ax1.plot(yvals,
             stats.beta.pdf(x=yvals, a=a, b=b),
             color=c[i], ls=ls[i], label=rf'$\alpha$={a}, $\beta$={b}')
    ax1.set(xlabel=r'確率変数 $\theta$', ylabel='density')
    ax1.legend(bbox_to_anchor=(0.5, -0.2), loc='upper center',
               title='パラメータ')
    ax1.grid(lw=0.3)
# 右のグラフの描画
for i, (a, b) in enumerate(zip([0.5, 0.5, 3.0], [0.5, 1.0, 1.0])):
    ax2.plot(yvals,
             stats.beta.pdf(x=yvals, a=a, b=b),
             color=c[i], ls=ls[i], label=rf'$\alpha$={a}, $\beta$={b:.1f}')
    ax2.set(xlabel=r'確率変数 $\theta$', ylabel='density', ylim=(-0.2, 5.2))
    ax2.legend(bbox_to_anchor=(0.5, -0.2), loc='upper center',
               title='パラメータ')
    ax2.grid(lw=0.3)
# 修飾
fig.suptitle('ベータ分布')
plt.tight_layout();

【実行結果】

6.5 カテゴリカル分布

図6.5を描画します。
確率分布には scipy.stats の rv_discrete() を利用します。

### 図6.5の描画

## 設定
# パラメータ
theta = [0.1, 0.2, 0.25, 0.35, 0.1]
# 色
c = 'tab:blue'

## 描画用データの作成
# パラメータKの算出
K = list(range(1, len(theta)+1))
# 確率質量関数の算出 ※汎用的な離散確率変数クラスの利用
# https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.rv_discrete.html
categorical = stats.rv_discrete(values=(K, theta))
prob = categorical.pmf(K)

## 描画処理
# 描画領域の設定
plt.figure(figsize=(5, 3))
ax = plt.subplot()
# 棒グラフの描画
ax.bar(K, prob, color=c)
# 修飾
ax.set(title=f'カテゴリカル分布\n$K={K},\ p={theta}$',
       xlabel='確率変数 $y$', ylabel='probability',
       xticks=K, ylim=(-0.02, 0.37), yticks=(0, 0.1, 0.2, 0.3))
plt.grid(lw=0.5);

【実行結果】

6.6 多項分布

テキストに掲載のない図を描きます。
確率分布には scipy.stats の multinomial() を利用します。

### 描画

## 設定
# パラメータ
N = 5
theta = [0.55, 0.45]  # K=2
# 色
c = 'tab:blue'

## 描画用データの作成
# yの算出
yvals = np.array([[i, N-i] for i in range(0, N+1)])
yvals_str = [str(yvals[i]) for i in range(len(yvals))]
# 確率質量関数の算出
prob = stats.multinomial.pmf(x=yvals, n=N, p=theta)

## 描画処理
# 描画領域の設定
plt.figure(figsize=(5, 3))
ax = plt.subplot()
# 棒グラフの描画
ax.bar(yvals_str, prob, color=c)
# 修飾
ax.set(title=f'多項分布\n$K={len(theta)},\ N={N},\ p={theta}$',
       xlabel='確率変数 $y$', ylabel='probability',
       ylim=(-0.02, 0.37), yticks=(0, 0.1, 0.2, 0.3))
plt.grid(lw=0.5);

【実行結果】

6.7 ディリクレ分布

図6.6を描画します。
確率分布には scipy.stats の dirichlet() を利用します。

三角図を自力で描画することができませんでした。
WolfMoon さんの関数コードをお借りしました。
ありがとうございます！

関数コードの定義です。後続処理に必要な箇所を改造しました。

### 三角図を描画する関数 by WolfMoon さんのコードを引用し、一部改造
# https://qiita.com/WolfMoon/items/0c90ec9f7645129859b7

def ternary_diagram(df, title='Ternary Diagram', lw=0.25, color='blue',
                    marker='o', s=10, alpha=0.5, index=False):
    """
    3 変数 `var1`、`var2`、`var3` の3列からなるデータフレームで与える
    （変数名はデータフレームの列名を転用する）。
    このプログラムは、正三角形の頂点座標を [0, 0],[1, 0],[1/2, np.sqrt(3) / 2] に
    設定し、`var1` を AB 方向の長さ、`var2` を BC 方向の長さ、`var3` を CA 方向の
    長さとして計算しプロットする。
    プロットに使う色は `color`，マーカーの種類は `marker`，サイズは `s` で
    指定できる。
    `index=True` を指定するとデータポイントの上に df.index をラベルとして表示する。
    使用例
    df = pd.DataFrame({'x' : [0.2, 0.3, 0.4],
                   'y' : [0.3, 0.2, 0.3],
                   'z' : [0.5, 0.5, 0.3]})
    ternary_diagram(df, title='三角グラフ')
    """

    # 正三角形の頂点座標（外枠を閉じるために4点目を追加）
    vertices = np.array([[0, 0], [1, 0], [1/2, np.sqrt(3) / 2], [0, 0]])

    # プロットデータの初期化
    n = len(df)
    labels = df.columns
    indexes = df.index
    plot_data = np.zeros((n, 2))
    var1 = df.iloc[:, 0].copy()
    var2 = df.iloc[:, 1].copy()
    var3 = df.iloc[:, 2].copy()        

    # 各変数の値を正規化してデータポイントの x-y 座標値を計算）
    for i in range(n):
        total = var1[i] + var2[i] + var3[i]
        var1[i] /= total
        var2[i] /= total
        var3[i] /= total
        plot_data[i] = (vertices[0] + var1[i] * (vertices[1] - vertices[0])
                        + var2[i] * (vertices[2] - vertices[0]))
        if index == True:
            ax.text(plot_data[i, 0], plot_data[i][1]+0.02, indexes[i],
                    ha='center', va='bottom', fontsize=8)


    # グリッド線の描画（AB, BC, CA に平行な破線）
    for a in [0.0, 0.2, 0.4, 0.6, 0.8, 1.0]:
        ax.text(1 - a/2, a * np.sqrt(3) / 2, f' {a}', ha='left', fontsize=9)
        ax.text(a/2, a * np.sqrt(3) / 2, f'{round(1 - a, 1)} ', ha='right',
                fontsize=9)
        ax.text(a, -0.05, a, ha='center', va='center', fontsize=9)
        if a != 0.0 and a != 1.0:
            ax.plot([a/2, 1 - a/2], [a * np.sqrt(3) / 2, a * np.sqrt(3) / 2],
                    'k--', linewidth=lw)
            ax.plot([a, a/2], [0, a * np.sqrt(3) / 2], 'k--', linewidth=lw)
            ax.plot([a, (1 + a)/2], [0, (1 - a) * np.sqrt(3) / 2], 'k--',
                    linewidth=lw)

    # プロット（外側の正三角形を描画，プロットポイントをマーカーで描画
    ax.plot(vertices[:, 0], vertices[:, 1], 'k-', linewidth=lw*2)
    ax.scatter(plot_data[:, 0], plot_data[:, 1], color=color, marker=marker,
               s=s, alpha=alpha)

    # 軸のラベルを追加
    delta = 0.125
    ax.text(0.5, -0.1, labels[0], ha='center', va='center', fontsize=12)
    ax.text(0.75+delta, np.sqrt(3)/4, labels[1], ha='center', va='center',
            fontsize=12, rotation=-60)
    ax.text(0.25-delta, np.sqrt(3)/4, labels[2], ha='center', va='center',
            fontsize=12, rotation=60)
    
    ax.set_xlim(-0.05, 1.05)
    ax.set_ylim(-0.15, 1.0)
    ax.set_title(title, fontsize=14)
    ax.axis('off')  # 座標軸を非表示
    ax.axis('equal')  # アスペクト比を1:1に設定

続いてディリクレ分布の確率分布を取得して、三角図を描画します。

### 図6.6の描画

## 設定：ディリクレ分布のパラメータαの値を設定
alphas = [[1, 1, 1], [0.3, 1, 1], [0.3, 0.3, 0.3], [3, 6, 6]]

## 描画処理
# 描画領域の指定
plt.figure(figsize=(8, 8))
# ４つのディリクレ分布のグラフを繰り返し描画
for i, alpha in enumerate(alphas):
    # ディリクレ分布の乱数を取得してデータフレーム化
    df = pd.DataFrame(
            stats.dirichlet.rvs(alpha=alpha, size=200, random_state=1),
            columns=['$\overrightarrow{θ}_1$', '$\overrightarrow{θ}_2$',
                     '$\overrightarrow{θ}_3$'])
    # 描画領域の指定
    ax = plt.subplot(2, 2, i+1)
    # 三角図の描画処理（関数利用）
    ternary_diagram(df, title=rf"$\overrightarrow{{θ}}$={alpha}", lw=0.25,
                    color='tab:blue', marker='o', s=20, alpha=0.5, index=False)
# 全体修飾
plt.suptitle('ディリクレ分布', fontsize=20)
plt.tight_layout();

【実行結果】

6.8 指数分布

図6.7を描画します。
確率分布には scipy.stats の expon() を利用します。

### 図6.7の描画

## 設定
# パラメータβ
betas = [1, 3]
# 色
c = ['tab:blue', 'tab:red']
# 線種
ls = ['-', '--']

## 描画用データの作成
# 確率変数yの値
yvals = np.linspace(0, 5, 1001)

## 描画処理
# 描画領域の設定
plt.figure(figsize=(5, 3))
ax = plt.subplot()
# 折れ線グラフの描画
for i, beta in enumerate(betas):
    ax.plot(yvals,
            stats.expon.pdf(x=yvals, scale=1/beta),
            color=c[i], ls=ls[i], label=rf'$\beta$={beta}')
# 修飾
ax.set(title='指数分布', xlabel='確率変数 $y$', ylabel='density')
ax.legend(bbox_to_anchor=(1.25, 1), title='パラメータ')
ax.grid(lw=0.5);

【実行結果】

6.9 ポアソン分布

図6.8を描画します。
確率分布には scipy.stats の poisson() を利用します。

### 図6.8の描画

## 設定
# パラメータλ
lam = 2.5
# 色
c = 'tab:blue'

## 描画用データの作成
# yの算出
yvals = np.arange(0, 11)
# 確率質量関数の算出
prob = stats.poisson.pmf(k=yvals, mu=lam)

## 描画処理
# 描画領域の設定
plt.figure(figsize=(5, 3))
ax = plt.subplot()
# 棒グラフの描画
ax.bar(yvals, prob, color=c)
# 修飾a
ax.set(title=f'ポアソン分布\n$\lambda={lam}$',
       xlabel='確率変数 $y$', ylabel='probability',
       xticks=range(0, 11), yticks=(0, 0.1, 0.2))
plt.grid(lw=0.5);

【実行結果】

6.10 ガンマ分布

図6.9を描画します。
確率分布には scipy.stats の gamma() を利用します。

### 図6.9の描画

## 設定
# 色
c = ['tab:blue', 'tab:red', 'tab:green']
# 線種
ls = ['-', '--', '-.']

## 描画用データの作成
# 確率変数yの値
yvals = np.linspace(0, 6, 1001)

## 描画処理
# 描画領域の設定
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(9, 5))
# 左のグラフの描画
for i, (a, b) in enumerate(zip([1, 3, 3], [1, 3, 1])):
    ax1.plot(yvals,
             stats.gamma.pdf(x=yvals, a=a, scale=1/b),
             color=c[i], ls=ls[i], label=rf'$\alpha$={a}, $\beta$={b}')
    ax1.set(xlabel='確率変数 $y$', ylabel='density',
            ylim=(-0.2, 2.1), yticks=np.arange(0, 2.1, 0.5))
    ax1.legend(bbox_to_anchor=(0.5, -0.2), loc='upper center',
               title='パラメータ')
    ax1.grid(lw=0.5)
# 右のグラフの描画
for i, (a, b) in enumerate(zip([1.0, 1.0, 0.1], [1.0, 3.0, 0.1])):
    ax2.plot(yvals,
             stats.gamma.pdf(x=yvals, a=a, scale=1/b),
             color=c[i], ls=ls[i], label=rf'$\alpha$={a}, $\beta$={b:.1f}')
    ax2.set(xlabel='確率変数 $y$', ylabel='density',
            ylim=(-0.2, 3.1), yticks=(0, 1, 2, 3),
            title=r'Gamma($\alpha$=1, $\beta$) = Exponential($\beta$)')
    ax2.legend(bbox_to_anchor=(0.5, -0.2), loc='upper center',
               title='パラメータ')
    ax2.grid(lw=0.5)
# 修飾
fig.suptitle('ガンマ分布')
plt.tight_layout();

【実行結果】

6.11 正規分布

図6.10を描画します。
確率分布には scipy.stats の norm() を利用します。

### 図6.10の描画

## 設定
# パラメータ μ, σ
mus = [0, 2]
sigmas = [1.0, 0.5]
# 色
c = ['tab:blue', 'tab:red']
# 線種
ls = ['-', '--']

## 描画用データの作成
# 確率変数yの値
yvals = np.linspace(-4, 4, 1001)

## 描画処理
# 描画領域の設定
plt.figure(figsize=(5, 3))
ax = plt.subplot()
# 折れ線グラフの描画
for i, (mu, sigma) in enumerate(zip(mus, sigmas)):
    ax.plot(yvals,
            stats.norm.pdf(x=yvals, loc=mu, scale=sigma),
            color=c[i], ls=ls[i], label=rf'$\mu$={mu}, $\sigma$={sigma}')
# 修飾
ax.set(title='正規分布', xlabel='確率変数 $y$', ylabel='density')
ax.legend(bbox_to_anchor=(1, 1), title='パラメータ')
ax.grid(lw=0.5);

【実行結果】

6.12 対数正規分布

図6.11を描画します。
確率分布には scipy.stats の lognorm() を利用します。

### 図6.11の描画

## 設定
# パラメータ μ, σ
mus = [0, 0, 2]
sigmas = [1.0, 0.4, 0.4]
# 色
c = ['tab:blue', 'tab:red', 'tab:green']
# 線種
ls = ['-', '--', '-.']

## 描画用データの作成
# 確率変数yの値
yvals = np.linspace(0, 10, 1001)

## 描画処理
# 描画領域の設定
plt.figure(figsize=(5, 3))
ax = plt.subplot()
# 折れ線グラフの描画
for i, (mu, sigma) in enumerate(zip(mus, sigmas)):
    ax.plot(yvals,
            stats.lognorm.pdf(x=yvals, s=sigma, scale=np.exp(mu)),
            color=c[i], ls=ls[i], label=rf'$\mu$={mu}, $\sigma$={sigma}')
# 修飾
ax.set(title='対数正規分布', xlabel='確率変数 $y$', ylabel='density')
ax.legend(bbox_to_anchor=(1.38, 1), title='パラメータ')
ax.grid(lw=0.5);

【実行結果】

6.13 多変量正規分布

図6.12を描画します。
確率分布には scipy.stats の multivariate_normal() を利用します。

### 図6.12の描画

## 描画用データ作成関数の定義
def make_data(mu_a, mu_b, sigma_a, sigma_b, rho, b_val):

       ## ２変量正規分布用のデータ生成
       # ２変量正規分布に用いるパラメータをくくる
       mean = [mu_a, mu_b]
       cov = [[sigma_a**2, rho * sigma_a * sigma_b],
              [rho * sigma_a * sigma_b, sigma_b**2]]
       # ２変量正規分布の乱数を200個生成
       random_num = stats.multivariate_normal.rvs(
                            mean=mean, cov=cov, size=200, random_state=0)
       
       ## 周辺分布 p(a)用のデータ生成
       xvals = np.linspace(-4, 6, 1001)
       yvals1 = stats.norm.pdf(xvals, loc=mu_a, scale=sigma_a)
       
       ## 条件付き分布 p(a | b=4)用のデータ生成 ※fig6-12.Rの計算式を引用
       mu_a_b = mu_a + rho * sigma_a / sigma_b * (b_val - mu_b)
       sigma_a_b = sigma_a * np.sqrt(1 - rho**2)
       yvals2 = stats.norm.pdf(xvals, loc=mu_a_b, scale=sigma_a_b)
       
       ## 戻り値
       # ２変量正規分布乱数、x軸の値、周辺分布確率密度、条件付き分布確率密度
       return random_num, xvals, yvals1, yvals2

## 左のグラフ用のデータ作成
# パラメータの設定 mu_a,b, sigma_a,b, rho
mu_a, sigma_a = 0, 1.5
mu_b, sigma_b = 1, 1.5
rho = 0.4
b_val = 4
# データ作成関数の実行
random_numL, xvalsL, yvals1L, yvals2L = make_data(
                            mu_a, mu_b, sigma_a, sigma_b, rho, b_val)

## 右のグラフ用のデータ作成
# パラメータの設定 mu_a,b, sigma_a,b, rho
mu_a, sigma_a = 1, 1.5
mu_b, sigma_b = 3, 0.5
rho = -0.7
b_val = 4
# データ作成関数の実行
random_numR, xvalsR, yvals1R, yvals2R = make_data(
                            mu_a, mu_b, sigma_a, sigma_b, rho, b_val)

### 描画処理
## 描画領域の設定
fig = plt.figure(figsize=(9, 5))
## ４つのグリッドの設定
gs = gridspec.GridSpec(2, 2, height_ratios=(1, 3))
ax1 = fig.add_subplot(gs[0, 0])
ax2 = fig.add_subplot(gs[1, 0])
ax3 = fig.add_subplot(gs[0, 1])
ax4 = fig.add_subplot(gs[1, 1])

## グリッド１つ目：左の周辺分布・条件付き分布の描画
ax1.plot(xvalsL, yvals1L)                            # 青：周辺分布
ax1.plot(xvalsL, yvals2L, color='tab:red', ls='--')  # 赤：条件付き分布
ax1.sharey(ax3)
ax1.grid(lw=0.5)

## グリッド２つ目：左の２変量正規分布乱数の散布図の描画
ax2.scatter(random_numL[:, 0], random_numL[:, 1], alpha=0.4)
ax2.axhline(4, color='tomato', ls='--', lw=2)
ax2.set(xlabel='a', ylabel='b')
ax2.sharex(ax1)
ax2.grid(lw=0.5)

## グリッド３つ目：右の周辺分布・条件付き分布の描画
ax3.plot(xvalsR, yvals1R, label='周辺分布')
ax3.plot(xvalsR, yvals2R, color='tab:red', ls='--', label='条件付き分布')
ax3.legend(bbox_to_anchor=(1, 1))
ax3.grid(lw=0.5)

## グリッド４つ目：右の２変量正規分布乱数の散布図の描画
ax4.scatter(random_numR[:, 0], random_numR[:, 1], alpha=0.4,
            label='２変量\n正規分布\n乱数')
ax4.axhline(4, color='tomato', ls='--', lw=2)
ax4.set(xlabel='a', ylabel='b')
ax4.sharex(ax3)
ax4.sharey(ax2)
ax4.legend(bbox_to_anchor=(1, 1))
ax4.grid(lw=0.5)

fig.suptitle('多変量正規分布')
plt.tight_layout()

【実行結果】

6.14 コーシー分布

図6.13を描画します。
確率分布には scipy.stats の cauchy() を利用します。

### 図6.13の描画

## 設定
# パラメータ μ, σ
mus = [0, 0]
sigmas = [1, 3]
# 色
c = ['tab:blue', 'tab:red']
# 線種
ls = ['-', '--']

## 描画用データの作成
# 確率変数yの値
yvals = np.linspace(-6, 6, 1001)

## 描画処理
# 描画領域の設定
plt.figure(figsize=(5, 3))
ax = plt.subplot()
# 折れ線グラフの描画
for i, (mu, sigma) in enumerate(zip(mus, sigmas)):
    ax.plot(yvals,
            stats.cauchy.pdf(x=yvals, loc=mu, scale=sigma),
            color=c[i], ls=ls[i], label=rf'$\mu$={mu}, $\sigma$={sigma}')
# 修飾
ax.set(title='コーシー分布', xlabel='確率変数 $y$', ylabel='density')
ax.legend(bbox_to_anchor=(1.35, 1), title='パラメータ')
ax.grid(lw=0.5);

【実行結果】

6.15 Student の t 分布

図6.14を描画します。
確率分布には scipy.stats の t() を利用します。

### 図6.14の描画

## 設定
# パラメータ ν, μ, σ
nus = [1, 4, np.inf]
mus = [0, 0, 0]
sigmas = [1, 1, 1]
# 色
c = ['tab:blue', 'tab:red', 'tab:green']
# 線種
ls = ['-', '--', '-.']

## 描画用データの作成
# 確率変数yの値
yvals = np.linspace(-6, 6, 1001)

## 描画処理
# 描画領域の設定
plt.figure(figsize=(5, 3))
ax = plt.subplot()
# 折れ線グラフの描画
for i, (nu, mu, sigma) in enumerate(zip(nus, mus, sigmas)):
    ax.plot(yvals,
            stats.t.pdf(x=yvals, df=nu, loc=mu, scale=sigma),
            color=c[i], ls=ls[i],
            label=rf'$\nu$={nu}, $\mu$={mu}, $\sigma$={sigma}')
# 修飾
ax.set(title='Studentのt分布', xlabel='確率変数 $y$', ylabel='density')
ax.legend(bbox_to_anchor=(1.45, 1), title='パラメータ')
ax.grid(lw=0.5);

【実行結果】

6.16 二重指数分布（ラプラス分布）

図6.15を描画します。
確率分布には scipy.stats の laplace() を利用します。

### 図6.15の描画

## 設定
# パラメータ μ, σ
mus = [0, 0]
sigmas = [1, 3]
# 色
c = ['tab:blue', 'tab:red']
# 線種
ls = ['-', '--']

## 描画用データの作成
# 確率変数yの値
yvals = np.linspace(-6, 6, 1001)

## 描画処理
# 描画領域の設定
plt.figure(figsize=(5, 3))
ax = plt.subplot()
# 折れ線グラフの描画
for i, (mu, sigma) in enumerate(zip(mus, sigmas)):
    ax.plot(yvals,
            stats.laplace.pdf(x=yvals, loc=mu, scale=sigma),
            color=c[i], ls=ls[i],
            label=rf'$\mu$={mu}, $\sigma$={sigma}')
# 修飾
ax.set(title='二重指数分布（ラプラス分布）', xlabel='確率変数 $y$',
       ylabel='density')
ax.legend(bbox_to_anchor=(1.45, 1), title='パラメータ')
ax.grid(lw=0.5);

【実行結果】

第６章の本論は以上です。

シリーズの記事

次の記事

前の記事

目次

ブログの紹介

note で７つのシリーズ記事を書いています。
ぜひ覗いていってくださいね！

１．のんびり統計

統計検定２級の問題集を手がかりにして、確率・統計をざっくり掘り下げるブログです。
雑談感覚で大丈夫です。ぜひ覗いていってくださいね。
統計検定２級公式問題集CBT対応版に対応しています。
Python、EXCELのサンプルコードの配布もあります。

２．実験！たのしいベイズモデリング１＆２をPyMC Ver.5で

書籍「たのしいベイズモデリング」・「たのしいベイズモデリング２」の心理学研究に用いられたベイズモデルを PyMC Ver.5で描いて分析します。
この書籍をはじめ、多くのベイズモデルは Ｒ言語＋Stanで書かれています。
PyMCの可能性を探り出し、手軽にベイズモデリングを実践できるように努めます。
身近なテーマ、イメージしやすいテーマですので、ぜひぜひPyMCで動かして、一緒に楽しみましょう！

３．実験！岩波データサイエンス1のベイズモデリングをPyMC Ver.5で

書籍「実験！岩波データサイエンスvol.1」の４人のベイジアンによるベイズモデルを PyMC Ver.5で描いて分析します。
この書籍はベイズプログラミングのイロハをざっくりと学ぶことができる良書です。
楽しくPyMCモデルを動かして、ベイズと仲良しになれた気がします。
みなさんもぜひぜひPyMCで動かして、一緒に遊んで学びましょう！

４．楽しい写経 ベイズ・Python等

ベイズ、Python、その他の「書籍の写経活動」の成果をブログにします。
主にPythonへの翻訳に取り組んでいます。
写経に取り組むお仲間さんのサンプルコードになれば幸いです🍀

５．RとStanではじめる心理学のための時系列分析入門 を PythonとPyMC Ver.5 で

書籍「RとStanではじめる心理学のための時系列分析入門」の時系列分析をPythonとPyMC Ver.5 で実践します。
この書籍には時系列分析のテーマが盛りだくさん！
時系列分析の懐の深さを実感いたしました。
大好きなPythonで楽しく時系列分析を学びます。

６．データサイエンスっぽいことを綴る

統計、データ分析、AI、機械学習、Pythonのコラムを不定期に綴っています。
統計・データサイエンス書籍にまつわる記事が多いです。
「統計」「Python」「数学とPython」「R」のシリーズが生まれています。

７．Python機械学習プログラミング実践記

書籍「Python機械学習プログラミング PyTorch & scikit-learn編」を学んだときのさまざまな思いを記事にしました。
この書籍は、scikit-learnとPyTorchの教科書です。
よかったらぜひ、お試しくださいませ。

最後までお読みいただきまして、ありがとうございました。

この記事が参加している募集

#仕事について話そう

112,434件

この記事が気に入ったらサポートをしてみませんか？