「データ解析のための統計モデリング入門」をPythonで写経 Vol.12 ～ ６章「GLMの応用範囲を広げる」④正規分布と尤度と分散パラメータ

６章「GLMの応用範囲を広げる」

書籍の著者　久保拓弥 先生

書籍「データ解析のための統計モデリング入門」６章「GLMの応用範囲を広げる」の Python写経活動記録 です。 

この記事は 正規分布の尤度 をきっかけにして、学びを広げていく様子を描いています。

テキスト６章が正規分布のGLMの統計モデリングを実践していない関係で、記事本編ではモデリングを行いません。
でもでも、アディショナルタイムで賑やかに寄り道します！
・正規分布のGLMモデリング

そして、ChatGPTの知識を借りてGLMモデリングの学びを深める過程で、素敵なことが分かったので、アディショナルタイムで共有させて下さい！
・GLMの分散パラメータの推定

では書籍を開いて統計モデリングの旅に出かけましょう🚀

はじめに

このブログシリーズは、書籍「データ解析のための統計モデリング入門 一般化線形モデル・階層ベイズモデル・MCMC」（岩波書店、「テキスト」と呼びます）の Python 写経を通じて得た「統計モデリングの楽しさ」をご紹介します。

テキストの紹介と引用表記はリンク先の記事に掲載しています。

準備

■ 記事の範囲
この記事はテキスト６章の以下の節を取り扱います。

■ ライブラリのインポート
この記事で用いるライブラリをインポートします。

# インポート

# 数値計算
import numpy as np

# 統計計算
import scipy.stats as stats

# 可視化
import matplotlib.pyplot as plt
plt.rcParams['font.family'] = 'Meiryo'  # または import japanize_matplotlib

正規分布と尤度

正規分布は確率変数が連続値をとる連続型確率分布です。

前回記事までに頻出したポアソン分布や二項分布は確率変数が離散値をとる離散型確率分布でした。
実は、連続型確率分布の尤度関数は離散型確率分布の尤度関数と意味合いが異なります。
テキストに沿って違いを確認していきます。

正規分布のプロフィール

正規分布の概要です。

【注記】パラメータ $${\sigma}$$ と $${\sigma^2}$$
テキストは正規分布の数式表現の際、ばらつきのパラメータに標準偏差 $${\sigma}$$ を用いています。
この記事では分散 $${\sigma^2}$$ を用います。

正規分布の確率密度関数

■ 数式表現
平均パラメータ $${\mu}$$、分散パラメータ $${\sigma^2}$$ の正規分布の確率密度関数は次の式で表されます。

「確率密度」は確率変数がある１点の値をとるときの「密度」であり、「確率」ではありません。

■ Python の関数表現
正規分布の確率密度関数を Python で実装してみます。
Python のライブラリの多くは、分散パラメータの代わりに標準偏差パラメータを用いています。
この記事も Python を使うときには標準偏差パラメータを使います。

# 正規分布の確率密度関数の定義
def norm_pdf(y, mu, sigma):
    return \
        1 / np.sqrt(2*np.pi * sigma**2) * np.exp(-(y - mu)**2 / (2 * sigma**2))

【実行結果】なし

確率変数 $${y=2}$$、平均パラメータ $${\mu=0}$$、標準偏差パラメータ $${\sigma=1}$$ として、正規分布の確率密度を計算します。

# 確率密度の計算 自作関数
y, mu, sigma = 2, 0, 1
norm_pdf(y=y, mu=mu, sigma=sigma)

【実行結果】
$${0.0540}$$ です。

scipy.stats の正規分布クラス norm の pdf メソッドで確率密度を計算します。

# 確率密度の計算 scipy.stats
y, mu, sigma = 2, 0, 1
stats.norm.pdf(x=y, loc=mu, scale=sigma)

【実行結果】
自作関数と同じ結果になりました。

■ 確率密度関数の描画
平均パラメータ $${\mu=0}$$、分散パラメータ $${\sigma^2=1}$$（標準偏差パラメータ $${\sigma=1}$$）の正規分布を標準正規分布と呼びます。
標準正規分布の確率密度関数の曲線を可視化します。

# 確率密度関数で標準正規分布の確率密度関数を可視化

## 設定
# 平均パラメータμ、標準偏差パラメータσ
mu, sigma = 0, 1
# x軸の確率変数yの値
y_vals = np.linspace(-4, 4, 1001)

## 描画
# 確率密度関数の描画
plt.plot(y_vals, norm_pdf(y_vals, mu, sigma), label='確率密度関数')
# y=2のときの確率密度の点と点線の描画
plt.plot(2, norm_pdf(2, mu, sigma), 'o', color='tab:red')
plt.vlines(2, -0.01, norm_pdf(2, mu, sigma), color='tab:red', ls='--')
# 修飾
plt.xlabel('確率変数 $y$', fontsize=12)
plt.ylabel('確率密度', fontsize=12)
plt.ylim(-0.01, None)
plt.show();

【実行結果】
正規分布の確率密度関数はベル型の形状です。
確率変数 $${y=2}$$ のときの確率密度に赤い点を描画しています。

正規分布の確率

■ 数式表現
連続型確率分布の「確率」は積分で定義されます。
正規分布に従う確率変数 $${y}$$が $${a}$$ 以上 $${b}$$ 以下となる確率 $${P(a \leq y \leq b)}$$ は、確率密度関数 $${f(y \mid \mu, \sigma^2)}$$ を用いて、次のように表されます。

■ 正規分布の確率の可視化
テキスト p.135 図 6.11 をお借りして、正規分布の「確率」を確認します。

# 正規分布の確率密度関数 p.135 図6.11

## 設定
mus, sigmas = [0, 0, 2], [1, 3, 1]        # 正規分布のパラメータ：μ, σ
y_vals = np.linspace(-5, 5, 200)          # x軸の値
y_vals_fill = np.linspace(1.2, 1.8, 100)  # x軸の値：塗りつぶし用
headers = ['A', 'B', 'C']                 # タイトルのかっこ内の文字

## 描画
# 描画領域の設定
fig, axes = plt.subplots(1, 3, figsize=(10, 3), tight_layout=True)
# ３つのグラフごとに正規分布の確率密度関数、青い領域の描画を繰り返し処理
for mu, sigma, header, ax in zip(mus, sigmas, headers, axes.flat):
    # 正規分布インスタンスの作成
    norm_dist_p = stats.norm(loc=mu, scale=sigma)
    # 正規分布の確率密度関数の描画: stats.norm.pdfで正規分布の確率密度関数を取得
    ax.plot(y_vals, norm_dist_p.pdf(y_vals))
    # 1.2≦y≦1.8の確率の塗りつぶし描画
    ax.fill_between(y_vals_fill, 0, norm_dist_p.pdf(y_vals_fill), alpha=0.3)
    # 確率の計算
    prob = norm_dist_p.cdf(1.8) - norm_dist_p.cdf(1.2)
    # 凡例
    ax.set(ylim=(0, 0.42), xlabel='y',
           title=f'({header}) $\mu$={mu}, $\sigma$={sigma}\n確率: {prob:.4f}')
axes[0].set_ylabel('確率密度');

【実行結果】
青い塗りつぶしの面積が確率です。
区間 $${[1.2, 1.8]}$$ の確率 $${\int_{1.2}^{1.8} f(y \mid \mu, \sigma^2) dy}$$ を表します。

尤度関数・対数尤度関数

■ 離散型確率分布の尤度関数
ポアソン分布等の離散型確率分布の尤度関数は、標本サイズ $${N}$$、パラメータ $${\theta}$$、確率質量関数 $${P(y \mid \theta)}$$ として、次のように「確率の積」で表されます。

テキストは、連続型確率分布である正規分布も「確率」で尤度関数を表すために、次のような工夫を用いて検討します。

■ 幅 $${\Delta y}$$ の導入
正規分布の確率を「図 (A) の面積計算をイメージして」計算します。
まず、確率変数が区間 $${[1.2, 1.8]}$$ の場合の標準正規分布の確率を計算します。

# 確率の計算 pnorm p.136

# μ=0, σ=1の正規分布のオブジェクトの定義
norm_dist = stats.norm(loc=0, scale=1)

# 1.2≦y≦1.8の確率の算出: cdfで累積分布関数を取得
norm_dist.cdf(x=1.8) - norm_dist.cdf(x=1.2)

【実行結果】

塗りつぶし領域を長方形に近似して面積を計算します。
確率変数が区間 $${[1.2, 1.8]}$$ の中間 $${1.5}$$ のときの確率密度を高さに、区間 $${[1.2, 1.8]}$$ の幅 $${\Delta y = 0.6}$$ にして、面積を計算します。

# 確率の近似計算 dnorm p.136
# 確率密度p(y=1.5|0,1)を高さ、1.8-1.2=0.6を幅Δyとする長方形で近似

norm_dist.pdf(x=1.5) * (1.8 - 1.2)

【実行結果】

テキストは幅 $${\Delta y}$$ が小さいほど良い近似になる、とコメントしています。

■  確率を使って正規分布の尤度関数・対数尤度関数を検討
確率 $${P(y_i - 0.5 \Delta y \leq y \leq y_i + 0.5 \Delta y)}$$ を確率密度関数と区間幅の積 $${f(y_i \mid \mu, \sigma^2) \Delta y}$$ で近似して、正規分布の尤度関数 $${L(\mu, \sigma^2)}$$ を書きます。

対数尤度関数は次のようになります。

◆ ◆ ◆

■ 正規分布の尤度関数・対数尤度関数
ここから世界が一変して、「連続型確率分布の尤度関数」の世界に飛びます！
定数である区間幅 $${\Delta y}$$ は、最尤推定におけるパラメータ $${\mu, \sigma}$$ の推定値に影響を与えないので、尤度関数・対数尤度関数から区間幅 $${\Delta y}$$ を外します。
この除外の結果、尤度関数と対数尤度関数は次のように変身します。

📊 正規分布の尤度関数・対数尤度関数

「確率」で尤度関数・対数尤度関数を検討をスタートしましたが、最終的に尤度関数・対数尤度関数は「確率密度」で構成されることになりました。

◆ ◆ ◆

■ 離散型確率分布と連続型確率分布の尤度関数の違い
離散型確率分布と連続型確率分布の尤度関数・対数尤度関数は「似ているけれども同じでは無い」ということが分かりました。
違いを表にまとめます。

またテキストは連続型確率分布の尤度の注意点を示しています。

正規分布のパラメータの最尤推定値

（注意）
こちらは趣味的な深堀りとコードです。
ご興味ない方はスルーしてくださって大丈夫です。

正規分布の対数尤度関数からパラメータ $${\mu, \sigma^2}$$ を推定します。
ここからは ChatGPT に全面依拠します！

以下では、正規分布の最尤推定について、「対数尤度 → 偏微分 → 最尤推定量の導出」をステップ・バイ・ステップでやさしく詳しく解説します。

🔧 前提：正規分布の確率密度関数
$${Y_i \sim \text{Normal}(\mu, \sigma^2)}$$ に従うときの確率密度関数は：

🧮 Step 1：尤度関数

独立な $${N}$$ 個の観測値 $${y_1, \dots, y_n}$$ に対する尤度関数 $${L(\mu, \sigma^2)}$$ は、各確率密度の積：

🧾 Step 2：対数尤度関数

対数を取って扱いやすく：

ここから $${\mu}$$ と $${\sigma^2}$$ に関して偏微分していきます。

✏️ Step 3：$${\mu}$$ について微分

これを 0 と置いて解くと：

➡️最尤推定値： $${\hat{\mu} = \bar{y}}$$（標本平均）

✏️ Step 4：$${\sigma^2}$$ について微分

これを 0 と置いて解くと：

➡️最尤推定値：$${\hat{\sigma}^2 = \frac{1}{N} \sum_i (y_i - \bar{y})^2}$$　（$${N}$$ で割る！）

✅ 結論：正規分布の最尤推定量

🔍 「$${N-1}$$」で割る不偏分散等ではなく、「$${N}$$」で割る点に注意。
最尤推定は分散の過小評価になることがありますが、これはバイアスとトレードオフで統計学的にはよく知られた結果です。

ChatGPT さん、ありがとうございました！
途中の計算に偏微分が入って難しかったですが、最尤推定値はとてもシンプルになりました。

ここまでは、分散 $${\sigma^2}$$ を最尤推定する前提のお話でした。
実はstatsmodels の GLM における分散は別の方法で推定されます。
最尤推定の「分散の過小評価」を回避できているのでご安心ください！
詳しくはアディショナルタイムで！

まとめ

今回は正規分布の尤度関数・対数尤度関数を学びました。

🔷 正規分布の尤度関数

🔷 正規分布の対数尤度関数

🔷 正規分布のパラメータの最尤推定値

アディショナルタイム

（注意）
アディショナルタイムの２テーマは趣味的な深堀りとコードです。
ご興味ない方はスルーしてくださって大丈夫です。

正規分布のGLMモデリング

統計モデルが無いと寂しいので（$${\approx}$$ ただやりたいので）、アディショナルタイムでパパパッとやってみます！
次回のガンマ分布の GLM で利用する例題データをお借りします。

■ 追加インポート

# 追加インポート
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf
import rdata

【実行結果】なし

■ データの読み込み
d.RData ファイルを rdata ライブラリで変換して、pandas データフレームの data に読み込みます。

# ガンマ分布のGLMで使用する例題データ

# RDataファイルの読み込み
data = rdata.read_rda('./data/ch06/d.RData', default_encoding='ASCII')['d']

# データフレームの表示
print('data.shape: ', data.shape)
data.head()

【実行結果】
データの個数（標本サイズ）は 50 です。
植物の個体（個体数 50）に関する仮想の観測データです。

【変数の説明】
植物の個体ごとの花重量 y と葉重量 x です。

◆ ◆ ◆

■ 確率分布、リンク関数、線形予測子

🔷 確率分布と確率密度関数
個体 $${i}$$ における花重量 $${y_i}$$ は平均花重量 $${\mu_i}$$、分散$${\sigma^2}$$ の正規分布に従います。

🔷 リンク関数と線形予測子
リンク関数は「恒等」、線形予測子は「$${\beta_1 + \beta_2 x_i}$$」です。

◆ ◆ ◆

■ Python ライブラリで統計モデルをデータに当てはめ
statsmodels の glm を利用して統計モデルを実装します。
引数 family で正規分布 Gaussian() を指定します。
リンク関数はデフォルトの恒等関数が適用されます。
当てはめ結果を変数 res_glm に格納します。

# 正規分布・恒等リンク関数のGLM（線形回帰）の当てはめ

# 設定
family = sm.families.Gaussian()

# モデルの当てはめ
res_glm = smf.glm(formula='y ~ x', family=family, data=data).fit()
res_glm.summary()

【実行結果】
最下２行の Intercept が $${\beta_1}$$、x が $${\beta_2}$$に対応しています。

切片 Intercept、葉重量 x の係数の推定値を線形予測子に当てはめてみます。

■ 最小二乗法 ols でフィッティング
正規分布・恒等リンク関数の GLM は線形回帰ですので、最小二乗法 OLS でモデルを推定することもできます。
statsmodels の ols でフィッティングします。
当てはめ結果を変数 res_ols に格納します。

# 最小二乗法olsでパラメータ推定すると…　パラメータ推定値の95%信頼区間が相違する
res_ols = smf.ols(formula='y ~ x', data=data).fit()
res_ols.summary()

【実行結果】

💡 statsmodels の GLM と OLS の相違点💡
パラメータ推定値とその標準誤差は同じ結果になりました。
しかし、統計量・$${p}$$ 値・95% 信頼区間が異なる値になっています。
理由は、GLM の統計量が標準正規分布に従う $${z}$$ 値、OLS の統計量が $${t}$$ 分布に従う $${t}$$ 値になっているからです。

■ OLS と GLM の使い分け（例）
ChatGPT に OLS と GLM の使い分けを訊いてみました。

✅ 一言アドバイス

• 「OLSとGLM（正規分布×恒等リンク）」はほぼ同じ係数を返しますが、周辺の扱い（標準誤差、信頼区間、推定の安定性）が異なるため、使い分けの意図をはっきり持つことが大事です。

• 教育目的や簡易分析には OLS、本格的な統計モデリングや将来的な拡張性が欲しいなら GLM が良い選択です。

線形回帰モデルの場合は、OLS を使うのがよさそうですね！

◆ ◆ ◆

■ 予測：最小二乗法 OLS の予測区間を使ってみる
OLS の結果を用いて、予測値の平均と 95% 信頼区間、新しい観測値に対する 95% 予測区間を可視化します。

# 平均μの予測値と95%信頼区間の可視化

## 設定と準備
# x軸の値
x_vals = np.linspace(0, 0.9, 101)
# olsの結果を用いて予測値のサマリーを取得
pred_ols = res_ols.get_prediction(dict(x=x_vals)).summary_frame()

## 描画
# x,yの散布図の描画
plt.plot(data.x, data.y, 'o', alpha=0.7)
# 平均μの予測値の描画
plt.plot(x_vals, pred_ols['mean'], color='tab:red', label='平均予測値（回帰直線）')
# 予測値の95%信頼区間の塗りつぶし描画
plt.fill_between(x_vals, pred_ols.mean_ci_lower, pred_ols.mean_ci_upper,
                 color='lightpink', alpha=0.5, label='平均予測の95%信頼区間');
# 新しい観測値に対する95%予測区間の塗りつぶし描画
plt.fill_between(x_vals, pred_ols.obs_ci_lower, pred_ols.obs_ci_upper,
                 color='lightpink', alpha=0.2, label='新しい観測値の95%予測区間');
# 修飾
plt.xlabel('葉重量 $x$', fontsize=12)
plt.ylabel('花重量 $y$', fontsize=12)
plt.legend();

【実行結果】
正規分布には下限値が無いため、信頼区間・予測区間にマイナス値が現れています。
著者の先生には、確率分布の仮定が大問題ですよ、と叱られそうです汗
また、薄赤色の 95% 信頼区間に含まれない観測値が２つあります。

【コード補足説明：予測区間の取得】

pred_ols = res_ols.get_prediction(dict(x=x_vals)).summary_frame()

最小二乗法の結果 res_ols に対して、get_prediction メソッドで予測を行い、さらに summary_frame メソッドで予測値関連情報を pandas データフレーム化しています。

summary_frame の結果を格納した pred_ols の中身を確認します。

# 予測値のサマリー: mean...:平均予測、obs...:新しい観測値に対する予測
pred_ols.head()

【実行結果】

平均予測値（$${\mu_i}$$）mean から始まり、平均予測の標準誤差 mean_se、平均予測の信頼区間 mean_ci_*、新しい観測値に対する予測区間 obs_ci_* で構成されています。

◆ ◆ ◆

■ GLM の結果でも予測区間を描画したい！
GLM の結果で予測区間を描画したいですか？👂️✋️
はい、とても描画したいです！
そこで、ChatGPT に正規分布の GLM の予測区間の算式を教えてもらいました。
$${t}$$ 分布を用いています。

📊 正規分布・恒等リンク関数の GLM の予測区間の公式

【変数の意味合い】

• $${x_0}$$：新たな説明変数（ベクトル、定数項を含む）

• $${y_0}$$：$${x_0}$$ に対応する新しい観測値（ベクトル）

• $${X}$$：モデリングに使用した説明変数のデザイン行列（定数項を含む）

• $${n}$$：標本サイズ（データ数）

• $${p}$$：パラメータ数（定数項を含む）

• $${t_{n-p,1-\alpha/2}}$$：自由度 $${n - p}$$ の t分布の上側 $${\alpha/2}$$ 点

では、Python 実装に移ります。
予測区間算出関数を定義します！

# 正規分布・恒等リンク関数のGLM用の予測区間算出関数

def pred_interval_for_gaussian_glm(result, x0, alpha=0.05):
    '''正規分布・恒等リンク関数のGLM用の予測区間算出関数
    
    Args:
        result (GLMResults): glmのresult（modelも包含していること）
        x0 (np.ndarray): 予測用の説明変数（先頭列に定数項を置く）
        alpha (float): 1 - 信頼係数（デフォルト 0.05）

    Returns:
        np.ndarray: 予測区間（1列目がlower, 2列目がupper）
    '''
    
    ## 設定と準備
    # 説明変数をnunpy配列化
    x0 = np.array(x0)
    # モデリングに使用したデザイン行列の取得
    X_design = result.model.exog
    # パラメータ数
    n_params = X_design.shape[1]
    
    ## 予測区間の算出
    # 計算式：
    # y_hat ± tval * se_pred
    # se_pred = sqrt(sigma2_hat * (1 + x0.T @ (X.T @ X).inv @ x0))

    # 平均予測値の算出
    y_hat = np.dot(x0, result.params)
    # 誤差の推定分散
    sigma2_hat = result.scale
    # (X.T @ X)^-1
    XtX_inv = np.linalg.inv(X_design.T @ X_design)
    # 標準誤差（予測区間用）
    se_pred = np.array([np.sqrt(sigma2_hat * (1 + x @ XtX_inv @ x)) for x in x0])
    # t値（自由度 = n - p）
    df = result.model.nobs - n_params
    tval = stats.t.ppf(1 - alpha / 2, df)
    # 予測区間の算出
    lower = y_hat - tval * se_pred
    upper = y_hat + tval * se_pred

    return np.column_stack([lower, upper])

【実行結果】なし

予測区間等の可視化に進みます！

# glmの結果を用いて、平均μの予測値と95%信頼区間の可視化

## 設定と準備
# x軸の値
x_vals = np.linspace(0, 0.9, 101)
# glmの結果を用いて予測値のサマリーを取得
pred_glm = res_glm.get_prediction(dict(x=x_vals)).summary_frame()
pred_int = pred_interval_for_gaussian_glm(
    res_glm, np.column_stack([np.ones_like(x_vals), x_vals]))
## 描画
# x,yの散布図の描画
plt.plot(data.x, data.y, 'o', alpha=0.7)
# 平均μの予測値の描画
plt.plot(x_vals, pred_glm['mean'], color='tab:red', label='平均予測値（回帰直線）')
# 予測値の95%信頼区間の塗りつぶし描画
plt.fill_between(x_vals, pred_glm.mean_ci_lower, pred_glm.mean_ci_upper,
                 color='lightpink', alpha=0.5, label='平均予測の95%信頼区間');
# 新しい観測値に対する95%予測区間の塗りつぶし描画
plt.fill_between(x_vals, pred_int[:, 0], pred_int[:, 1],
                 color='lightpink', alpha=0.2, label='新しい観測値の95%予測区間');

# 修飾
plt.xlabel('葉重量 $x$', fontsize=12)
plt.ylabel('花重量 $y$', fontsize=12)
plt.legend();

【実行結果】
OLS の結果を用いた予測区間と同じような帯になりました。

このGLM の予測区間が OLS の予測区間と一致しているか確かめます。
両者の予測区間の下側値が 1e-16 の精度で一致していることを確認するコードです。

# olsの予測区間（下側）とglmの予測区間（下側）の計算結果の一致
np.allclose(pred_ols.obs_ci_lower, pred_int[:, 0], rtol=0, atol=1e-16)

【実行結果】
1e-16 の精度で一致しました（1e-17 だと一部不一致になります）。

GLMの分散パラメータの推定

ところで、正規分布と言えば平均パラメータ $${\mu}$$ と分散パラメータ $${\sigma^2}$$ が特徴的な確率分布ですよね。
GLM のパラメータ推定では、$${\mu}$$ に関連する係数 $${\beta_1, \beta_2}$$ の推定値を中心に見てきました。
そういえば、$${\sigma^2}$$ の推定値を見掛けていませんね…

正規分布のGLM の結果を再掲します。

下部のパラメータ推定値のゾーンに分散は見当たりません！

でもご安心ください。
実は、上部の「Scale」が分散の推定値なのです！
この分散パラメータ（dispersion parameter）を $${\phi}$$ としましょう。

ChatGPT によると…

なんと！
GLM の裏側（からくり）が少し分かった気がします。

$${\phi}$$ の計算方法を ChatGPT と Gemini に教わりましょう。

🎯 正規分布の GLM における分散パラメータ $${\bm \phi}$$ 
正規分布の GLM の場合、$${\phi}$$ は観測値 $${y_i}$$ の分散 $${\text{Var}(y_i)}$$ です。

◆

🧮 $${\bm \phi}$$ の推定式 
statsmodels における GLM の分散パラメータ $${\phi}$$ は、正規分布などの一部の確率分布の場合、デフォルトで、ピアソンカイ二乗統計量（Pearson's chi-squared statistic）を残差の自由度で割った値として推定されます。
この値は、$${\mathtt{fit()}}$$ メソッド実行後に $${\mathtt{result.scale}}$$ として確認できます。

📐 ピアソンカイ二乗統計量に基づく算出方法（一般論）

• $${\hat{\phi}}$$：分散パラメータの推定値

• $${N}$$：観測数

• $${p}$$：モデルのパラメータ数

• $${y_i}$$：$${i}$$番目の観測値

• $${\hat{\mu}_i}$$：予測された平均値

• $${V(\hat{\mu}_i)}$$：分散関数（確率分布に依存）

📐 ピアソンカイ二乗統計量に基づく算出方法（正規分布のGLMの場合）
正規分布の場合、$${V(\hat{\mu}_i)=1}$$ で計算します。

GLM の深淵を覗いた気がします…（感動）
正規分布以外を含めて、GLM 全般の分散パラメータ $${\phi}$$ を訊きました。

📊 $${\mathtt{scale}}$$（$${\phi}$$）のデフォルト計算式の詳細

◆

⚙️ $${\mathtt{fit()}}$$ メソッドの $${\mathtt{scale}}$$ 引数 
statsmodels の GLM は設定によって、分散パラメータの計算方法を変えることができます。
具体的には、$${\mathtt{fit()}}$$ の $${\mathtt{scale}}$$ 引数で指定します。
何も指定しないとデフォルトが適用されます。

• $${\mathtt{scale=None}}$$（デフォルト）：正規分布・ガンマ分布・逆正規分布の場合は $${\mathtt{scale='X2'}}$$、二項分布・ポアソン分布の場合は $${\mathtt{scale=1.0}}$$ を適用

• $${\mathtt{scale='X2'}}$$：ピアソンカイ二乗統計量を残差の自由度で割って算出

• $${\mathtt{scale=1.0}}$$：ポアソンや二項分布など、分散パラメータを1に固定したい場合

• 任意の数値を指定：特定の値に固定したいとき

• $${\mathtt{scale='dev'}}$$：逸脱度を残差の自由度で割って算出

◆

🧭 注意点 
ポアソン分布や二項分布のように、分散が平均に依存する「分散パラメータが1に固定される」ファミリーでは、通常 $${\mathtt{scale}}$$ は1として扱われます。
ただし、データに過分散（overdispersion）がある場合、statsmodels は自動的に $${\mathtt{scale}}$$ を1より大きい値に推定し、その影響を考慮することができます。
これは、$${\mathtt{fit()}}$$ メソッドで $${\mathtt{scale='X2'}}$$ を指定した場合などに起こり、よりロバストな標準誤差や $${p}$$ 値が得られます。

ChatGPT & Gemini ありがとうございます！

先ほどの正規分布の GLM の結果から $${\mathtt{scale}}$$ 属性で分散パラメータを取り出してみましょう。

# 分散パラメータの推定値の表示
res_glm.scale

【実行結果】
$${\hat{\phi} = 0.0107}$$ です。

ピアソンカイ二乗統計量を残差の自由度で割って分散パラメータを計算します。
ピアソンカイ二乗統計量は $${\mathtt{result.pearson\_chi2}}$$、残差の自由度は $${\mathtt{result.df\_resid}}$$ で取得できます。

# 誤差の推定分散の計算
# 誤差の推定分散 = ピアソンカイ二乗統計量 / 残差の自由度
res_glm.pearson_chi2 / res_glm.df_resid

【実行結果】
分散パラメータ $${\phi}$$ の値になりました！

実は正規分布の GLM の場合、別の方法で分散パラメータを計算することができます。
引き続き、ChatGPT と Gemini に教わりましょう。

🧮 $${\bm \phi}$$ の推定式（正規分布の場合の特別な計算方法） 
正規分布を使ったGLMでは、分散パラメータ（$${\phi}$$）は、残差平方和（RSS）を残差の自由度で割った値と等しくなります。

◆

🧠 背景
正規分布のGLMでは、分散関数 $${V(\mu)}$$ は定数 1 になります。
これは、正規分布の分散が平均 $${\mu}$$に依存しないためです。

• 正規分布の分散：$${\mathrm{Var}(y) = \sigma^2}$$

• GLMとしての分散：$${\mathrm{Var}(y) = \phi \cdot V(\mu) = \phi \cdot 1 = \phi}$$

したがって、正規分布の GLM では、分散パラメータ $${\phi}$$ がそのまま正規分布の分散（誤差分散） $${\sigma^2}$$ を表します。

◆

📐 ピアソンカイ二乗統計量からの導出
GLMにおける分散パラメータの推定式：

正規分布の場合、$${V(\hat{\mu}_i) = 1}$$ のため：

◆

✅ 結論
この結果は、正規分布・恒等リンク関数によるGLMが、通常の最小二乗法（OLS）と等価であることと一致します。
つまり、OLS は GLMの一部として自然に位置づけられ、$${\mathtt{scale}}$$ は OLS の分散推定量と一致します。

重回帰モデルなどの「最小二乗法」と正規分布・恒等リンク関数の GLM が結びつきました！

この残差平方和を用いた分散パラメータの計算をやってみます。
残差は $${\mathtt{result.resid\_response}}$$、残差の自由度は $${\mathtt{result.df\_resid}}$$ で取得できます。

# 誤差の推定分散の計算
# 誤差の推定分散 = 残差平方和 / 残差の自由度
sum(res_glm.resid_response**2) / res_glm.df_resid

【実行結果】
$${\hat{\phi} = 0.0107}$$ です。
先ほどの計算値と一致しています！

今回のブログは以上です。

次回はガンマ分布の GLM を学びます。

シリーズの記事

次の記事

前の記事

目次

ブログの紹介

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の教科書です。
よかったらぜひ、お試しくださいませ。

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