「データ解析のための統計モデリング入門」をPythonで写経 Vol.13 ～ ６章「GLMの応用範囲を広げる」⑤ガンマ分布と対数リンク関数のGLM
６章「GLMの応用範囲を広げる」
書籍の著者 久保拓弥 先生
書籍「データ解析のための統計モデリング入門」６章「GLMの応用範囲を広げる」の Python写経活動記録 です。
この記事は ガンマ分布と対数リンク関数のGLM を学びます。
ガンマ分布は正規分布と同じ「連続型確率分布」です。
では書籍を開いて統計モデリングの旅に出かけましょう🚀
はじめに
このブログシリーズは、書籍「データ解析のための統計モデリング入門 一般化線形モデル・階層ベイズモデル・MCMC」（岩波書店、「テキスト」と呼びます）の Python 写経を通じて得た「統計モデリングの楽しさ」をご紹介します。
テキストの紹介と引用表記はリンク先の記事に掲載しています。
準備
準備
■ 記事の範囲
この記事はテキスト６章の以下の節を取り扱います。
6.8 ガンマ分布の GLM
■ 利用データ
テキスト・サポートサイトのデータファイルを引用しています。
▶️ サポートサイト
Jupyter Notebook ファイルと同一フォルダ内に「data」フォルダを用意して、data フォルダ配下の章別フォルダにデータファイルを格納しています。
■ ライブラリのインポート
この記事で用いるライブラリをインポートします。
# インポート
# 数値計算
import numpy as np
import pandas as pd
# 統計計算
import scipy.stats as stats
import statsmodels.api as sm
import statsmodels.formula.api as smf
# Rデータセットの読み込み
import rdata
# 可視化
import matplotlib.pyplot as plt
import seaborn as sns
plt.rcParams['font.family'] = 'Meiryo' # または import japanize_matplotlib
統計モデリング・サマリー
この記事で扱う統計モデリングの概要です。
■ 統計モデル
ガンマ分布と対数リンク関数を用いる統計モデルです。
$$
\begin{array}{clll}
確率分布 & リンク関数 & モデル名 & パラメータ数 k \\
\hline
\\
ガンマ分布 & 対数 & \mathtt{x + f} モデル & k=3 \\
\end{array}
$$
■ モデリング手続き
1️⃣データの確認
2️⃣ガンマ分布
3️⃣統計モデルをデータに当てはめ
・統計モデルの理解
・当てはめと評価
4️⃣予測
ガンマ分布＋対数リンク関数のGLM
データの確認
データを読み込み、データの外観を眺めてから、統計モデリングのためのデータの特徴確認を行います。
■ RDataファイルの読み込み
d.RData ファイルを rdata ライブラリで変換して、pandas データフレームの data に読み込みます。
### ガンマ分布を使ったGLMの例題 p.139~140 図6.13
# RDataファイルの読み込み
data = rdata.read_rda('./data/ch06/d.RData', default_encoding='ASCII')['d']
# csvファイルの出力
# data.to_csv('./data/ch06/data_d.csv')
# データフレームの表示
print('data.shape: ', data.shape)
data.head()【実行結果】
データの個数（標本サイズ）は 50 です。
植物の個体（個体数 50）に関する仮想の観測データです。
【変数の説明】
植物の個体ごとの花重量 y（目的変数）と葉重量 x（説明変数）です。
$$
\begin{array}{clll}
変数 & 説明 & 値 \\
\hline
\\
y & 個体 i の花重量 & 0以上の実数 \\
x & 個体 i の葉重量 & 0以上の実数 \\
\end{array}
$$
◆ ◆ ◆
■ データの確認
基本的な統計量やチャートでデータを概観します。
① 要約統計量の表示
# 要約統計量
data.describe().round(3)【実行結果】
上からデータの個数、平均、標準偏差、最小値、第１四分位数、中央値、第三四分位数、最大値です。
標準偏差は不偏分散の標準偏差です。
② 標本分散の表示
要約統計量に含まれない分散を確認します。
# 標本分散
data.var(ddof=1).rename('var').to_frame().T.round(3)【実行結果】
花重量 y の標本分散は標本平均 0.177 よりも一桁小さな値です。
③ 相関係数の表示
# 相関係数
data.corr().round(3)【実行結果】
葉重量 x と花重量 y との間に中程度の正の相関があります。
⑤ ヒストグラムと散布図の描画
seaborn の pairplot() を利用します。
# ヒストグラムと散布図の描画
sns.pairplot(data=data, height=3,
diag_kws={'edgecolor': 'white', 'alpha': 0.7});【実行結果】
【考察】
花重量 y のヒストグラムは左に偏って右に裾の長い分布を示しています。実はガンマ分布は、右に裾の長い特徴を持っています。
葉重量 x のヒストグラムは一様分布を示しています。
横軸 x、縦軸 y の散布図（左下）は、葉重量 x が大きくなるにつれて、葉重量 y が大きくなる傾向があります。
しかも、y が大きくなるにつれて縦幅（ばらつき）も増大しています。
実はガンマ分布の分散は平均の関数になっています。
◆ ◆ ◆
■ データの特徴まとめ
データの特徴を整理します。
① 花重量は０以上の実数（連続値）
② 花重量の分布は右に裾が長い
③ 花重量は値が大きくなるにつれて、ばらつきが大きくなる
④ 葉重量と花重量は正の相関がある
データのばらつきは「確率分布」で表現します。
上述のデータの特徴を表現できる確率分布の候補は「ガンマ分布」です。
今回のデータは０以上の連続値、右に裾の長い分布、値が大きくなるにつれてばらつきが大きくなるなどの特徴がガンマ分布と合っていそうです。
ガンマ分布
ガンマ分布は確率変数が連続値をとる連続型確率分布です。
パラメータには複数の定義があるようです。
この記事では３つのパラメータ定義をご紹介いたします。
1️⃣ shape、rate パラメータ
１つ目の定義はテキストが採用する定義です。
R 言語で使われているようです。
📶 ガンマ分布の確率密度関数
確率変数$${y}$$、形状（shape）パラメータ $${s}$$、率（rate）パラメータ $${r}$$、ガンマ関数 $${\Gamma(\cdot)}$$ を用いて、ガンマ分布の確率密度関数は次の式で表されます。
$$
f(y \mid s, r) = \cfrac{r^s}{\Gamma(s)}\ y^{s-1} \exp(-ry) \\
$$
SciPy 公式サイトの $${\mathtt{gamma}}$$ クラスのページに掲載の
$$
f(x \mid \alpha, \beta) = \cfrac{\beta^{\alpha}}{\Gamma(\alpha)} x^{\alpha - 1} \exp(-\beta x)
$$
のように、形状パラメータ $${\alpha}$$、率パラメータ $${\beta}$$ と表記するケースもあります。
📶 ガンマ分布のプロフィール
$$
\begin{array}{lll}
項目 & 内容 \\
\hline
\\
パラメータ & 形状パラメータ\ s>0& \\
& 率パラメータ\ r > 0 \\
確率変数 & [0, \inf) の実数 \\
& ただし\text{GLM} では (0, \inf) \\
平均 & s / r & \\
分散 & s / r^2 & \\
\end{array}
$$
📶 ガンマ分布の確率密度関数の可視化
テキスト p.139 図 6.12 に相当します。
確率変数 $${y}$$ の範囲 $${[0, 5]}$$、パラメータ $${s=r=[1, 5, 0.1]}$$ のケースで描画します。
区間 $${[1.2, 1.8]}$$ の確率（面積）を併せて描画します。
# ガンマ分布の確率密度関数 p.139 図6.12
## 設定
r_ss = [0.1, 1, 5] # ガンマ分布分布のパラメータ：r, s
y_vals = np.linspace(0, 5, 200) # x軸の値
y_vals_fill = np.linspace(1.2, 1.8, 100) # x軸の値：塗りつぶし用
headers = ['C', 'A', 'B'] # タイトルのかっこ内の文字
## 描画
# 描画領域の設定
fig, axes = plt.subplots(1, 3, figsize=(10, 4), sharey=True, tight_layout=True)
# ３つのグラフごとにガンマ分布の確率密度関数、青い領域の描画を繰り返し処理
for r_s, header, ax in zip(r_ss, headers, axes.flat):
# ガンマ分布のパラメータの取得
r = s = r_s
# ガンマ分布の確率密度関数の描画
# stats.gamma.pdfでガンマ分布の確率密度関数を取得
# パラメータの読み替え:a=s, scale=1/r
gamma_dist = stats.gamma(a=s, scale=1 / r)
ax.plot(y_vals, gamma_dist.pdf(y_vals))
# 1.2≦y≦1.8の確率の塗りつぶし描画
ax.fill_between(y_vals_fill, 0, gamma_dist.pdf(y_vals_fill), alpha=0.3)
# x軸=0の垂直線の描画
ax.axvline(0, color='black', ls='--', lw=0.5)
# 修飾
prob = gamma_dist.cdf(1.8) - gamma_dist.cdf(1.2)
ax.set_title(f'({header}) $r=s=${r_s}\n確率={prob:.4f}')
ax.set(ylim=(0, 1.1), xlabel='y')
axes[0].set_ylabel('確率密度');【実行結果】
左に偏り、右に裾の長い分布です。
$${s, r}$$ の値が大きくなるにつれて峰が現れます。
$${s=1}$$ のとき、指数分布になります。
◆ ◆ ◆
2️⃣ shape、scale パラメータ
２つ目の定義は Scipy のパラメータが採用する定義です。
Wikipedia のガンマ分布もこの定義です。
📶 ガンマ分布の確率密度関数
確率変数$${y}$$、形状（shape）パラメータ $${\alpha}$$、尺度（scale）パラメータ $${\beta}$$、ガンマ関数 $${\Gamma(\cdot)}$$ を用いて、ガンマ分布の確率密度関数は次の式で表されます。
$$
f(y \mid \alpha, \beta) = \cfrac{1}{\Gamma(\alpha) \beta^{\alpha}}\ y^{\alpha-1} \exp(-y/\beta) \\
$$
📶 ガンマ分布のプロフィール
$$
\begin{array}{lll}
項目 & 内容 \\
\hline
\\
パラメータ & 形状パラメータ\ \alpha>0& \\
& 尺度パラメータ\ \beta > 0 \\
確率変数 & [0, \inf) の実数 \\
& ただし\text{GLM} では (0, \inf) \\
平均 & \alpha \beta & \\
分散 & \alpha \beta^2 & \\
\end{array}
$$
📶 １つ目の定義との関係
形状パラメータは同じです。
率パラメータと尺度パラメータは逆数の関係です。
$$
\begin{align*}
\alpha &= s \\
\beta &= 1/ r \\
r &= 1/ \beta \\
\end{align*}
$$
📶 ガンマ分布の確率密度関数の可視化
テキスト p.139 図 6.12 に相当します。
確率変数 $${y}$$ の範囲 $${[0, 5]}$$、パラメータ $${s=r=[1, 5, 0.1]}$$ のケースを $${\alpha, \beta}$$ に変換して描画します。
区間 $${[1.2, 1.8]}$$ の確率（面積）を併せて描画します。
# ガンマ分布の確率密度関数 p.139 図6.12 形状パラメータα, 率パラメータβ版
## 設定
r_ss = [0.1, 1, 5] # ガンマ分布分布のパラメータ：r, s
y_vals = np.linspace(0, 5, 200) # x軸の値
y_vals_fill = np.linspace(1.2, 1.8, 100) # x軸の値：塗りつぶし用
headers = ['C', 'A', 'B'] # タイトルのかっこ内の文字
## 描画
# 描画領域の設定
fig, axes = plt.subplots(1, 3, figsize=(10, 4), sharey=True, tight_layout=True)
# ３つのグラフごとにガンマ分布の確率密度関数、青い領域の描画を繰り返し処理
for r_s, header, ax in zip(r_ss, headers, axes.flat):
# ガンマ分布のパラメータα,βの取得
alpha = r_s
beta = 1 / r_s # β = 1 / r
# ガンマ分布の確率密度関数の描画
# stats.gamma.pdfでガンマ分布の確率密度関数を取得
# パラメータの読み替え:a=α, scale=β
gamma_dist = stats.gamma(a=alpha, scale=beta)
ax.plot(y_vals, gamma_dist.pdf(y_vals))
# 1.2≦y≦1.8の確率の塗りつぶし描画
ax.fill_between(y_vals_fill, 0, gamma_dist.pdf(y_vals_fill), alpha=0.3)
# x軸=0の垂直線の描画
ax.axvline(0, color='black', ls='--', lw=0.5)
# 修飾
prob = gamma_dist.cdf(1.8) - gamma_dist.cdf(1.2)
ax.set_title(f'({header}) $\\alpha=${alpha}, $\\beta=${beta}'
f'\n確率={prob:.4f}')
ax.set(ylim=(0, 1.1), xlabel='y')
axes[0].set_ylabel('確率密度');【実行結果】
左に偏り、右に裾の長い分布です。
$${\alpha, \beta}$$ の値と分布の関係は、イマイチ掴み取りにくいですね汗
ChatGPT によると…
$${\alpha < 1}$$：右に大きく歪んだ尖った分布（指数分布に近い）
$${\alpha = 1}$$：指数分布（exp. dist）
$${\alpha > 1}$$：山を持ち、右裾が長い分布
$${\alpha \gg 1}$$：正規分布に近づく（中心極限定理）
◆ ◆ ◆
3️⃣ 平均、ばらつきパラメータ
３つ目の定義は GLM と相性がよい定義です。
📶 ガンマ分布の確率密度関数
確率変数$${y}$$、平均パラメータ $${\mu}$$、分散パラメータ（dispersion parameter） $${\phi}$$、ガンマ関数 $${\Gamma(\cdot)}$$ を用いて、ガンマ分布の確率密度関数は次の式で表されます。
$$
f(y \mid \mu, \phi) = \cfrac{1}{\Gamma(1 / \phi)} \left( \cfrac{1}{\phi \mu} \right)^{1/\phi} \ y^{1 / \phi - 1} \exp \left(- \cfrac{y}{\phi\mu} \right) \\
$$
📶 ガンマ分布のプロフィール
$$
\begin{array}{lll}
項目 & 内容 \\
\hline
\\
パラメータ & 平均パラメータ\ \mu>0& \\
& 分散パラメータ\ \phi > 0 \\
確率変数 & [0, \inf) の実数 \\
& ただし\text{GLM} では (0, \inf) \\
平均 & \mu & \\
分散 & \phi \mu^2 & \\
\end{array}
$$
📶 １つ目・２つ目の定義との関係
$$
\begin{align*}
\mu &= \cfrac{s}{r} = \alpha \beta \\
\phi &= \cfrac{1}{s} = \cfrac{1}{\alpha} \\
\end{align*}
$$
📶 ガンマ分布の確率密度関数の可視化
テキスト p.139 図 6.12 に相当します。
確率変数 $${y}$$ の範囲 $${[0, 5]}$$、パラメータ $${s=r=[1, 5, 0.1]}$$ のケースを $${\mu, \phi}$$ に変換して描画します。
区間 $${[1.2, 1.8]}$$ の確率（面積）を併せて描画します。
# ガンマ分布の確率密度関数 p.139 図6.12 平均パラメータμ, 分散パラメータφ版
## 設定
r_ss = [0.1, 1, 5] # ガンマ分布分布のパラメータ：r, s
y_vals = np.linspace(0, 5, 200) # x軸の値
y_vals_fill = np.linspace(1.2, 1.8, 100) # x軸の値：塗りつぶし用
headers = ['C', 'A', 'B'] # タイトルのかっこ内の文字
## 描画
# 描画領域の設定
fig, axes = plt.subplots(1, 3, figsize=(10, 4), sharey=True, tight_layout=True)
# ３つのグラフごとにガンマ分布の確率密度関数、青い領域の描画を繰り返し処理
for r_s, header, ax in zip(r_ss, headers, axes.flat):
# ガンマ分布のパラメータμ,φの取得
mu = r_s / r_s # μ = s / r
phi = 1 / r_s # φ = 1 / s
# ガンマ分布の確率密度関数の描画
# stats.gamma.pdfでガンマ分布の確率密度関数を取得
# パラメータの読み替え:a=1/φ, scale=φμ
gamma_dist = stats.gamma(a=1 / phi, scale=phi * mu)
ax.plot(y_vals, gamma_dist.pdf(y_vals))
# 1.2≦y≦1.8の確率の塗りつぶし描画
ax.fill_between(y_vals_fill, 0, gamma_dist.pdf(y_vals_fill), alpha=0.3)
# x軸=0の垂直線の描画
ax.axvline(0, color='black', ls='--', lw=0.5)
# 修飾
prob = gamma_dist.cdf(1.8) - gamma_dist.cdf(1.2)
ax.set_title(f'({header}) $\mu=${mu}, $\phi=${phi}\n確率={prob:.4f}')
ax.set(ylim=(0, 1.1), xlabel='y')
axes[0].set_ylabel('確率密度');【実行結果】
左に偏り、右に裾の長い分布です。
$${\mu, \phi}$$ の値と分布の関係は、イマイチ掴み取りにくいですね汗
そこで ChatGPT の協力のもと、パラメータ値の変化を掴みやすそうな可視化を試みます。
# 平均パラメータμ, 分散パラメータφ とガンマ分布
## 設定
# 描画用のμ、φの組み合わせ
params = [
(1, 1, '基準'),
(2, 1, '右にシフト'),
(0.5, 1, '左にシフト'),
(1, 2, '広がる'),
(1, 0.1, '尖る'),
]
# x軸の値
x = np.linspace(1e-6, 4, 500)
## 描画
# 描画領域の設定
plt.figure(figsize=(8, 5))
# μ,φの組み合わせごとにガンマ分布の確率密度関数の描画を繰り返し処理
for mu, phi, label in params:
# scipy.stats.gamma用のパラメータ値の設定
shape = 1 / phi
scale = phi * mu
# ガンマ分布の確率密度関数の算出
y = stats.gamma.pdf(x, a=shape, scale=scale)
# 描画
plt.plot(x, y, label=f'$\mu$={mu}, $\phi$={phi}（{label}）')
# 修飾
plt.title('ガンマ分布 $\mu, \phi$ パターン比較')
plt.xlabel('$y$', fontsize=12)
plt.ylabel('確率密度', fontsize=12)
plt.ylim(0, 2.5)
plt.grid(alpha=0.5)
plt.legend()
plt.show()【実行結果】
✨ グラフの読み取りポイント
$${\mu}$$ を動かして $${\phi}$$ を固定
$${\mu}$$ が大きくなるにつれて、分布全体が右にスライド
（中心が右に）$${\phi}$$ を動かして $${\mu}$$（中心）を固定
$${\phi}$$ が小さいと分布は峰が高くなって尖る
$${\phi}$$ が大きくなるにつれて、分布は広がって平らになる
◆ ◆ ◆
■ Python でガンマ分布の確率密度関数を表現
① 確率密度関数の数式を関数化
３つの定義を包含して関数化します。
$${y=3, \alpha=2, \beta=3}$$ で関数の出力をテストします。
# ガンマ分布の確率密度関数の実装
from scipy.special import gamma as Gamma
# ガンマ分布の確率密度関数の定義 ※gamma
def gamma_pdf(y, param1, param2, kind=3):
match kind:
case 1: # 1. テキスト・Rと同じ定義 param1 = s, param2 = r
shape = param1
scale = 1 / param2
var_names = ('s', 'r')
case 2: # 2. scipyと同じ定義 param1 = α, param2 = β
shape = param1
scale = param2
var_names = ('α', 'β')
case 3: # 3. GLMと同じ定義 param1 = mu, param2 = phi
shape = 1 / param2
scale = param1 * param2
var_names = ('μ', 'φ')
# 戻り値：ガンマ分布の確率密度関数、パラメータ名
return (
(y**(shape - 1) * np.exp(- y / scale)) / (Gamma(shape) * scale**shape),
var_names
)
# テスト
# y, shape, scale, kind = 3, 2, 1/3, 1 # ケース1 R
y, shape, scale, kind = 3, 2, 3, 2 # ケース2 scipy
# y, shape, scale, kind = 3, 6, 1/2, 3 # ケース3 GLM
pdf, var_names = gamma_pdf(y, shape, scale, kind)
print(f'y={y}, {var_names[0]}={shape}, '
f'{var_names[1]}={scale} のガンマ分布確率密度 f(y) = {pdf:.4f}')【実行結果】
② scipy.stats の ガンマ分布クラス利用
gamma.pdf() で確率密度関数を算出できます。
# scipy で計算
y, shape, scale = 3, 2, 3
stats.gamma.pdf(x=y, a=shape, scale=scale)【実行結果】
◆ ◆ ◆
■ 観測データとガンマ分布の特徴
観測データの特徴はガンマ分布と似ています。
$$
\begin{array}{ll}
データの特徴 & ガンマ分布 \\
\hline
\\
0以上の実数（連続値） & 確率変数が0以上の実数 \\
\\
分布は右に裾が長い & パラメータによっては \\
&分布は右に裾が長い \\
\\
値が大きくなるにつれて & 分散は平均の二乗に比例 \\
ばらつきが大きくなる & （３つ目の定義）\\
\end{array}
$$
ガンマ分布を観測データに当てはめる際の「前提条件」をテキストから引用いたします。
・データが連続値
・データが $${[0, \infty]}$$
・分散は平均の関数
統計モデルをデータに当てはめ（モデルの理解）
今回はガンマ分布・対数リンク関数を用いる GLM を例題データ＝観測データに当てはめます。
モデル名は「ガンマ分布モデル」（非公式）です。
パラメータの定義は「３つ目」を使います。
◆ ◆ ◆
■ 確率分布、リンク関数、線形予測子
GLMの３要素である「確率分布」、「リンク関数」、「線形予測子」を定義します。
🔷 確率分布と確率密度関数
花重量 $${y_i}$$ は、平均パラメータ $${\mu_i}$$、分散パラメータ $${\phi}$$ のガンマ分布に従います。
$$
\begin{align*}
y_i &\sim \text{Gamma}(\mu_i, \phi) \\
f(y_i \mid \mu_i, \phi) &= \cfrac{1}{\Gamma(1 / \phi)} \left( \cfrac{1}{\phi \mu_i} \right)^{1/\phi} \ y^{1 / \phi - 1} \exp \left(- \cfrac{y}{\phi\mu_i} \right)
\end{align*}
$$
🔷 リンク関数と線形予測子
リンク関数は「対数」、線形予測子は「$${a + b \log x_i}$$」です。
$$
\begin{align*}
\log \mu_i &= a + b \log x_i \\
\Longleftrightarrow \mu_i &= \exp(a)x_i^b = \exp(a + b \log x_i)
\end{align*}
$$
◆
■ 平均 $${\mu_i}$$ を可視化して、形状のイメージアップへ
$${\mu_i = \exp(a)x_i^b = \exp(a + b \log x_i)}$$ を可視化します。
$${\mu_i = \exp(a)x_i^b}$$ に注目して、$${\exp(a)=1}$$（$${a=0}$$）、$${b=0.5}$$ の曲線を描きます。
$${x_i^{0.5} = \sqrt{x_i}}$$ です。
比較要素として、$${x_i^2}$$ を併せて描画します。
# x_i^b のグラフの描画
# x軸の値の設定
x_vals = np.linspace(0, 1, 100)
# x^0.5=√x の曲線の描画
plt.plot(x_vals, x_vals**0.5, color='tab:blue', label='$x^{0.5}=\sqrt{x}$')
# x^2 の曲線の描画
plt.plot(x_vals, x_vals**2, color='gray', ls='--', label='$x^2$')
# 修飾
plt.xlabel('$x_i$', fontsize=14)
plt.ylabel('$\mu_i$', fontsize=14)
plt.legend();【実行結果】
青い実線が $${x_i^{0.5} = \sqrt{x_i}}$$ の曲線です。
テキスト p.139 図 6.13 の曲線によく似ていることを確認できます。
$${x_i^b}$$ の指数 $${b}$$（上の例では $${0.5}$$） と、$${x_i^b}$$ に掛ける係数 $${\exp(a)}$$ の $${a}$$（上の例では $${0}$$）を推定します。
統計モデルをデータに当てはめ（当てはめと評価）
■ 統計モデルをデータに当てはめ、の準備
GLMの３要素「確率分布」「リンク関数」「線形予測子」を用いる統計モデルをデータに当てはめします。
統計モデルの対数尤度 $${\log L}$$ が最大になるパラメータ $${\beta_1, \beta_2, \beta_3}$$（線形予測子のパラメータ）を推定します。
🔷 今回の統計モデルの尤度関数 $${L}$$ と対数尤度関数 $${\log L}$$
$$
\begin{align*}
L(a, b) &= \prod_{i=1}^{n} \left(
\cfrac{1}{\Gamma(1/\phi)} \left( \cfrac{1}{\phi \mu_i} \right)^{1/\phi}
y_i^{1/\phi - 1} \exp\left( -\cfrac{y_i}{\phi \mu_i} \right)
\right)
\\
\\
\log L(a, b) &= \sum_{i=1}^{n} \left(
-\log\Gamma(1/\phi)
+ \cfrac{1}{\phi} \log\left( \frac{1}{\phi \mu_i} \right)
+ \left( \cfrac{1}{\phi} - 1 \right) \log y_i
- \cfrac{y_i}{\phi \mu_i}
\right) \\
\\
\mu_i &= \exp(a + b \log x_i) \\
\end{align*}
$$
GLM のパラメータ推定では、線形予測子のパラメータ（今回は $${a, b}$$）を最尤推定したあとに、別途、$${\phi}$$ を推定するそうです。
ですので、$${L, \log L}$$ の変数に $${\phi}$$ を含めていません。
なお、最適化問題（たとえば、scipy の minimize）で $${\phi}$$ も最尤推定する場合には、$${\phi}$$ を尤度関数・対数尤度関数の変数に含めることになります。
◆ ◆ ◆
■ Python ライブラリで統計モデルをデータに当てはめ
statsmodels の glm を利用して統計モデルを実装します。
引数 family にはガンマ分布 sm.families.Gamma() を設定します。
引数 link にはリンク関数 Log() を設定します。
当てはめ結果を変数 result に格納します。
### ガンマ分布を使ったGLMの例題 p.139~140 図6.13
# GLMモデルの当てはめ: 目的変数のばらつき⇒ガンマ分布、リンク関数⇒対数
## 設定
# formulaの設定
formula = 'y ~ np.log(x)'
# ガンマ分布・対数リンク関数の設定
family = sm.families.Gamma(link=sm.families.links.Log())
## モデリング
# モデルの当てはめ glm(y ~ log(x), family=Gamma(link='log'), data=d)
result = smf.glm(formula=formula, family=family, data=data).fit()
# 結果表示
result.summary()【実行結果】
最下２行の Intercept が $${a}$$、np.log(x) が $${b}$$ に対応しています。
◆ ◆ ◆
■ 当てはめ結果の分析
🔷 係数の推定値
切片 Intercept、葉重量 x の係数の推定値を線形予測子に当てはめてみます。
$$
\mu_i = \exp(-1.0403 + 0.6832 \log x_i)
$$
この数式を用いて、$${\mu}$$ の予測値を可視化しましょう。
# muの予測値グラフの描画
## 設定
# x軸の値の設定
x_vals = np.linspace(1e-6, 1, 100)
## muの予測値の算出
# 線形予測子のパラメータa,bの取得
a, b = result.params
# 平均muの予測値の計算
mus = np.exp(a + b * np.log(x_vals))
## 描画
# 観測値の散布図の描画
plt.plot(data.x, data.y, 'o', color='tab:blue', alpha=0.7)
# muの予測値の描画
plt.plot(x_vals, mus, color='tab:red', label='$\mu$ の予測値')
# x^0.5=√x の曲線の描画
plt.plot(x_vals, x_vals**0.5, color='gray', label='$x^{0.5}=\sqrt{x}$')
# 修飾
plt.xlabel('$x_i$', fontsize=14)
plt.ylabel('$\mu_i$', fontsize=14)
plt.legend();【実行結果】
赤い実線が $${\mu}$$ の予測値です。
さきほど描画した $${x_i^{0.5} = \sqrt{x_i}}$$（グレイ） よりも係数 $${a, b}$$ の値が小さいことが分かります。
🔷 パラメータ推定値の評価
２つ係数の推定値の $${p}$$ 値（P>|z|）は $${0.000<0.05}$$、95% 信頼区間（[0.025 0.975]）はゼロを含んでいません。
統計的な視点では、パラメータ推定値は有意ということになります。
また、係数の推定値の標準誤差は、推定値と比べて小さな値であり、特段の問題は無いように思われます。
🔷 モデルの評価
AIC で予測の良さを確認します。
# Null Deviance, Residual Deviance, AICの表示
print(f'Null 逸脱度\t: {result.null_deviance:5.1f}')
print(f'残差逸脱度\t: {result.deviance:5.1f}')
print(f'AIC\t\t: {result.aic:5.1f}')【実行結果】
AIC は -112.9 です。マイナス値になりました。
比較するモデルが無いので、ここでは参考値の扱いです。
◆ ◆ ◆
■ 最適化ライブラリで最尤推定
こちらは趣味のコードです。
ご興味ない方はスルーしてくださって大丈夫です。
statsmodels を使わないで最尤推定を行ってみます。
scipy の最小化問題ソルバー mininize で、負の対数尤度が最小となるパラメータを推定します。
# 最適化ライブラリで最尤推定
# 追加インポート
from scipy.optimize import minimize
# 目的関数の定義：負の対数尤度関数 ※負の対数尤度の最小化問題を解く
def llf(params, x, y):
# パラメータの分解 ※phiは対数にして正値を保つ
a, b, log_phi = params
phi = np.exp(log_phi)
# μ = exp(線形予測子)
mu = np.exp(a + b * np.log(x))
# ガンマ分布のパラメータの算出
shape = 1 / phi
scale = phi * mu
# 戻り値：負の対数尤度
return -np.sum(stats.gamma.logpdf(x=y, a=shape, scale=scale))
# 最適化の実行
res_optim = minimize(fun=llf, x0=[0, 0, 1], args=(data['x'], data['y']),
tol=1e-20)
# 結果の表示
print(f'a の最尤推定値 = {res_optim.x[0]:9.6f}')
print(f'b の最尤推定値 = {res_optim.x[1]:9.6f}')
print(f'φ の最尤推定値 = {np.exp(res_optim.x[2]):9.6f}')
print(f'最大対数尤度 = {-res_optim.fun:9.6f}')【実行結果】
パラメータ $${a, b}$$ は statsmodels の当てはめ結果と同じになっています。
一方で…というか、やはり、分散パラメータ $${\phi}$$ はズレました。
GLM の分散パラメータ推定値と minimize による分散パラメータの最尤推定値を比べます。
# 分散パラメータφ の比較
print(f'GLMのscale\t: {result.scale:.6f}')
print(f'最尤推定値\t: {np.exp(res_optim.x[2]):.6f}')【実行結果】
両者は異なる値になっています。
分散パラメータの推定方法が GLM と 最尤推定で異なるので、両者の分散パラメータ推定値が異なっても仕方がないのです。
◆ ◆ ◆
■ GLMの分散パラメータあれこれ
GLM は確率分布で目的変数のバラツキ（誤差）を表現しています。
このバラツキの素が分散パラメータ $${\phi}$$ です。
分散パラメータの推定値は、GLMの結果の scale 欄に表示されます。
今回は scale 欄に0.32508 と表示されています。
🧮 ガンマ分布の分散パラメータ $${\phi}$$ の推定式
$$
\hat{\phi} = \frac{1}{n - p} \sum_{i=1}^{n} \left( \frac{y_i - \hat{\mu}_i}{\hat{\mu}_i} \right)^2
$$
$${n}$$：データ数（サンプル数）
$${p}$$：推定したパラメータ数（回帰係数の数）
$${y_i}$$：観測値
$${\hat{\mu}_i}$$：予測された平均値（逆リンク関数を通したもの）
この推定式を使って、分散パラメータを算出してみます。
# ガンマ分布のGLMの分散パラメータ
# 計算要素の取得
n = result.model.nobs
p = len(result.params)
y = result.model.endog
mu_hat = result.mu
# ガンマ分布のGLMの分散パラメータφの算出
1 / (n - p) * sum(((y - mu_hat) / mu_hat)**2)【実行結果】
分散パラメータを算出できました。
GLM の分散パラメータは、ピアソンカイ二乗統計量を残差の自由度で割って求めます。
上の推定式の $${\sum_{i=1}^{n} \left( \frac{y_i - \hat{\mu}_i}{\hat{\mu}_i} \right)^2}$$ がピアソンカイ二乗統計量に相当します。
GLM の結果からピアソンカイ二乗統計量と残差の自由度を取得して確かめてみます。
# GLMの分散パラメータφの算出 ※ピアソンカイ二乗統計量 / 残差の自由度
result.pearson_chi2 / result.df_resid【実行結果】
GLM の結果に示された scale の値 0.32508と一致しています。
花重量の予測
テキストにならって、花重量 y の予測を行って可視化します。
■ 真のパラメータを既知として、予測・予測区間を可視化
GLM の結果 result に対して predict メソッドを適用して予測の平均を算出しています。
予測区間は「真のパラメータ値が分かっている前提で」サポートサイトのRスクリプトのロジックを模して算出しています。
p.139 図 6.13 に相当します。
### ガンマ分布を使ったGLMの例題 p.139~140 図6.13
## 設定と準備
# 設定
b1, b2 = -1, 0.7 # 真のパラメータ値
x_vals = np.linspace(data.x.min(), data.x.max(), 100) # 予測に用いるx軸の値
true_means = np.exp(b1 + b2 * np.log(x_vals)) # 真の平均
# 予測区間のためのパラメータ推定値の算出 ※テキストのRスクリプトを活用
phi = result.scale # 分散パラメータ
shape = 1 / phi # shapeパラメータ
rate = 1 / (phi * true_means) # rateパラメータ
## 描画
# 描画領域の設定
fig, ax = plt.subplots()
# 観測値の散布図の描画
sns.scatterplot(data=data, x='x', y='y', s=70, alpha=0.7, label='観測値',
zorder=10, ax=ax)
# 真の平均の曲線の描画
ax.plot(x_vals, true_means, color='black', ls='--', label='真の平均')
# 平均の予測値の曲線の描画
ax.plot(x_vals, result.predict(dict(x=x_vals)), color='royalblue',
label='予測の平均')
# 予測分布の中央値の曲線の描画
ax.plot(x_vals, stats.gamma.ppf(q=0.5, a=shape, scale=1 / rate), color='tab:red',
label='予測分布の中央値')
# 50%予測区間の塗りつぶし描画 ※y軸の区間の算出はテキストのRスクリプトを活用
ax.fill_between(x_vals,
stats.gamma.ppf(q=0.25, a=shape, scale=1 / rate),
stats.gamma.ppf(q=0.75, a=shape, scale=1 / rate),
color='lightpink', alpha=0.5, label='50%予測区間')
# 90%予測区間の塗りつぶし描画 ※y軸の区間の算出はテキストのRスクリプトを活用
ax.fill_between(x_vals,
stats.gamma.ppf(q=0.05, a=shape, scale=1 / rate),
stats.gamma.ppf(q=0.95, a=shape, scale=1 / rate),
color='lightpink', alpha=0.2, label='90%予測区間')
# 修飾
ax.set_xlabel('葉重量 $x_i$', fontsize=12)
ax.set_ylabel('花重量 $y_i$', fontsize=12)
plt.legend();【実行結果】
葉重量が大きくなるにつれて、花重量は大きくなり、かつ、花重量のバラツキも大きくなっている様子がよく分かります！
ちなみにテキスト p.140 脚注 *45 に「予測区間を評価するときに、推定値の誤差も考慮する場合があります。ここでは考慮していません」と補足されています。
この脚注について、ChatGPT は次のようにコメントしています。
書籍の脚注が示していたのはまさに「パラメータ推定値の誤差（不確かさ）を考慮していない」という点で、この予測区間は「 $${\phi}$$ を固定して、観測誤差のみを含めた予測分布による区間」の特徴そのものでした。
✅ では本題：推定値の誤差も考慮した予測区間とは？
これは、次の2つの不確かさを同時に考慮した予測区間です：
・パラメータ（$${\hat{a}, \hat{b}}$$）の推定誤差
・観測誤差（ガンマ分布のばらつき）
このような予測区間を構成するためには、以下のような方法があります：
✅ 方法①：ブートストラップ法
✅ 方法②：ベイズ的予測分布（posterior predictive）
✅ 方法③：近似的な修正（デルタ法＋ガンマ分布）
なるほど、ですね。
後ほどブートストラップ法を試してみます。
◆ ◆ ◆
■ 真のパラメータを未知として、予測・予測区間を可視化
その前に、真のパラメータを「未知」として、予測区間を可視化してみます。
$${\mathtt{rate}}$$ パラメータ（逆数は $${\beta}$$ パラメータ）の計算に用いる平均に、真の平均（$${\mathtt{true\_means}}$$）の代わりに、推定した平均（つまり予測の平均）を使います。
# ガンマ分布を使ったGLMの例題
# 真のパラメータ値が分からない場合、平均予測値で代用する
## 設定と準備
# 設定
x_vals = np.linspace(data.x.min(), data.x.max(), 100) # 予測に用いるx軸の値
x0 = dict(Intercept=np.ones_like(x_vals), x=x_vals) # 予測に用いる説明変数
mu_hat = result.predict(x0) # μの平均予測値
# 予測区間のためのパラメータ推定値の算出
phi = result.scale # 分散パラメータ
shape = 1 / phi # shapeパラメータ
scale = phi * mu_hat # scaleパラメータ
## 描画
# 描画領域の設定
fig, ax = plt.subplots()
# 観測値の散布図の描画
sns.scatterplot(data=data, x='x', y='y', s=70, alpha=0.7, label='観測値',
zorder=10, ax=ax)
# 平均予測値の曲線の描画
ax.plot(x_vals, result.predict(dict(x=x_vals)), color='royalblue',
label='予測の平均')
# 予測分布の中央値の曲線の描画
ax.plot(x_vals, stats.gamma.ppf(q=0.5, a=shape, scale=scale), color='tab:red',
label='予測分布の中央値')
# 50%予測区間の塗りつぶし描画
ax.fill_between(x_vals,
stats.gamma.ppf(q=0.25, a=shape, scale=scale),
stats.gamma.ppf(q=0.75, a=shape, scale=scale),
color='lightpink', alpha=0.5, label='50%予測区間')
# 90%予測区間の塗りつぶし描画
ax.fill_between(x_vals,
stats.gamma.ppf(q=0.05, a=shape, scale=scale),
stats.gamma.ppf(q=0.95, a=shape, scale=scale),
color='lightpink', alpha=0.2, label='90%予測区間')
# 修飾
ax.set_xlabel('葉重量 $x_i$', fontsize=12)
ax.set_ylabel('花重量 $y_i$', fontsize=12)
plt.legend();【実行結果】
真の平均を使った前のチャートと比べて、予測区間の幅が狭くなった気がします。
◆ ◆ ◆
こちらは趣味のコードです。
ご興味ない方はスルーしてくださって大丈夫です。
■ 真のパラメータを未知、予測区間にパラメータ推定誤差を含める
ブートストラップ法を用いて、パラメータ推定誤差と観測誤差の両方を含める予測区間を描画しましょう。
ブートストラップ法で「予測値のサンプルデータ」を生成する関数を定義します。
観測データから復元抽出法でランダムに標本サイズと同じ数の「ブートストラップ」サンプルを取得して、このサンプルにガンマ分布の GLM を当てはめ⇒推定したパラメータと $${\hat{\mu}}$$ を取得⇒ガンマ分布乱数（観測値の予測）の生成、を行います。
コードの骨格は ChatGPT によります。
# ガンマ分布・対数リンク関数のGLMでブートストラップ予測サンプル生成関数の定義
def gamma_glm_bootstrap_preds(data, formula, X0, B=1000, random_state=42):
'''
ガンマ分布・対数リンク関数のGLMの統計モデルを用いて、
複数の予測点 X0 に対して、y0のブートストラップ予測サンプルを一括生成する
Parameters:
data : pd.DataFrame of shape（N, P）
学習に用いる説明変数（列名付き）
formula : str
statsmodels用のGLMモデル式（例：'y ~ np.log(x)'）
X0 : dataframe-like of shape (N, P)
予測に用いる説明変数（列名付き）
B : int
ブートストラップ回数
random_state : int
乱数シード
Returns:
samples : ndarray of shape (N, B)
各予測点に対する B 個の予測値サンプル
'''
## インポート
import numpy as np
import statsmodels.formula.api as smf
from tqdm.notebook import tqdm # プログレスバー
## 設定と準備
# 乱数生成器
rng = np.random.default_rng(random_state)
# 予測値の個数
N = len(X0)
# 予測値のサンプルを格納する配列の初期化
samples = np.empty((N, B))
## ブートストラップ予測サンプルの生成
for b in tqdm(range(B)):
# ブートストラップリサンプリング
sample = data.sample(n=len(data), replace=True)
# ブートストラップサンプルデータにモデルを当てはめ
result = smf.glm(
formula,
data=sample,
family=sm.families.Gamma(sm.families.links.Log())
).fit()
# 予測平均値 μ_hat の算出（統計モデルで予測）
mu_hat = result.predict(X0) # shape (N,)
# ガンマ分布のパラメータの算出
phi_hat = result.scale
shape = 1 / phi_hat
scale = phi_hat * mu_hat # shape (N,)
# ガンマ分布から予測値をサンプリング
samples[:, b] = rng.gamma(shape=shape, scale=scale)
return samples【実行結果】なし
5000 個の予測値サンプルデータを生成します。
# ブートストラップ予測サンプルの生成
## 設定と準備
# 予測に用いるXの値
x_vals_bs = pd.Series(np.linspace(data.x.min(), data.x.max(), 21), name='x')
# GLMのformula
formula = result.model.formula
## ブートストラップ予測サンプル生成の実行
samples = gamma_glm_bootstrap_preds(data, formula, x_vals_bs, B=5000)【実行結果】
処理時間は 40 秒強でした。
では予測区間の描画に進みます！
# ガンマ分布を使ったGLMの例題 推定値の誤差を含む予測の可視化
## 設定と準備
# 予測サンプルの分位数を取得するヘルパー関数の定義
quantiles = lambda q: np.quantile(samples, q=q, axis=1)
# 平均予測値の算出
pred_gamma = result.predict(x_vals_bs)
## 描画
# 描画領域の設定
fig, ax = plt.subplots()
# 観測値の散布図の描画
sns.scatterplot(data=data, x='x', y='y', s=70, alpha=0.7, label='観測値',
zorder=10, ax=ax)
# 平均予測値の曲線の描画
ax.plot(x_vals_bs, pred_gamma, color='royalblue', label='予測の平均')
# 予測分布の中央値の曲線の描画
ax.plot(x_vals_bs, quantiles(q=0.5), color='tab:red', label='予測分布の中央値')
# 50%予測区間の塗りつぶし描画
ax.fill_between(x_vals_bs, *quantiles(q=[0.25, 0.75]),
color='lightpink', alpha=0.5, label='50%予測区間')
# 90%予測区間の塗りつぶし描画
ax.fill_between(x_vals_bs, *quantiles(q=[0.05, 0.95]),
color='lightpink', alpha=0.2, label='90%予測区間')
# 修飾
ax.set_xlabel('葉重量 $x_i$', fontsize=12)
ax.set_ylabel('花重量 $y_i$', fontsize=12)
plt.legend();【実行結果】
１つ前の「真のパラメータは未知」「予測区間にパラメータ推定誤差を含まない」予測区間とあまり変わらない結果になりました…
まとめ
今回はガンマ分布の GLM を実践しました。
パラメータの定義には複数あり、このまとめは、３つ目の GLM と相性の良いパラメータを振り返ります。
🔷 確率分布と確率密度関数
$$
\begin{align*}
y_i &\sim \text{Gamma}(\mu_i, \phi) \\
f(y_i \mid \mu_i, \phi) &= \cfrac{1}{\Gamma(1 / \phi)} \left( \cfrac{1}{\phi \mu_i} \right)^{1/\phi} \ y^{1 / \phi - 1} \exp \left(- \cfrac{y}{\phi\mu_i} \right)
\end{align*}
$$
🔷 リンク関数と線形予測子
$$
\begin{align*}
\log \mu_i &= a + b \log x_i \\
\Longleftrightarrow \mu_i &= \exp(a)x_i^b = \exp(a + b \log x_i)
\end{align*}
$$
🔷 statsmodels のロジスティック回帰モデル構築と結果表示
formula = 'y ~ np.log(x)'
family = sm.families.Gamma(link=sm.families.links.Log())
result = smf.glm(formula=formula, family=family, data=data).fit()
result.summary()今回のブログは以上です。
次回は、一般化線形混合モデル（GLMM）を実践します。
シリーズの記事
次の記事
前の記事
目次
ブログの紹介
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の教科書です。
よかったらぜひ、お試しくださいませ。
最後までお読みいただきまして、ありがとうございました。
いいなと思ったら応援しよう！
応援ありがとうございます。これからもがんばって記事を作成します！