「データ解析のための統計モデリング入門」をPythonで写経 Vol.2 ～ ２章「確率分布と統計モデルの最尤推定」 統計モデリング、起動

２章「確率分布と統計モデルの最尤推定」

書籍の著者　久保拓弥 先生

書籍「データ解析のための統計モデリング入門」２章「確率分布と統計モデルの最尤推定」の Python写経活動記録 です。 

この記事は シンプルな統計モデルを用いて統計モデリングを起動 します。

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

はじめに

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

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

統計モデリング、起動

準備

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

■ 利用データ
テキスト・サポートサイトのデータファイルを引用しています。
▶️ サポートサイト

Jupyter Notebook ファイルと同一フォルダ内に「data」フォルダを用意して、data フォルダ配下の章別フォルダにデータファイルを格納しています。

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

### インポート

# 数値計算
import math
import numpy as np
import pandas as pd

# 統計計算
import scipy.stats as stats

# Rデータセットの読み込み
import rdata

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

統計モデリング・サマリー

この記事で扱う統計モデリングの概要です。

■ 統計モデル

■ モデリング手続き

データの確認

データを読み込み、データの外観を眺めてから、統計モデリングのためのデータの特徴確認を行います。

■ RDataファイルの取り扱い
テキストは R 言語を利用しており、配布データは R 特有の RData 形式で作成されています。
rdata ライブラリで、RData 形式のデータを Python で扱えるデータ形式に変換します。

## データの読み込み

# RdataをPythonオブジェクト（今回は辞書型）に変換
converted = rdata.read_rda('./data/ch02/data.RData', default_encoding='ASCII')
converted

【実行結果】

【コード解説】
rdata ライブラリの read_rdata 関数でPython のデータ形式に変換します。

RData ファイルの文字コードが ASCII のようですので、引数 default_encoding で 'ASCII' を指定しました。
実行結果を見ると Python の辞書型（key='data'）に変換できたことが分かります。

■ pandas データフレーム化
辞書型のデータは pandas のデータフレームに簡単に取り込めます。
列名を「種子数」、データ型を整数型（int）に設定してデータフレーム data を作成します。

# pandasデータフレーム化
data = pd.DataFrame(converted['data'], columns=['種子数']).astype(int)

# 結果の表示
print('data.shape: ', data.shape)
display(data.head())

# 読み込んだデータをcsvファイルで保存
# data.to_csv('./data/ch02/data.csv', encoding='utf_8_sig', index=False)

【実行結果】
データの個数（標本サイズ）は 50、列は整数値の種子数データです。

◆ ◆ ◆

■ データの確認
テキスト p.15～17 のデータ要約・可視化を通じて、データの特徴を確認します。
テキストの初出しデータということで、データ操作に関する説明が厚い（熱い）です。

① data の内容の表示
種子数列を values 属性で numpy 配列に変換すると、テキストのように横並びに表示できます。

# dataの内容表示 p.15
data['種子数'].values 

【実行結果】

② データの個数（標本サイズ）の表示
len 関数でデータの個数をカウントできます。

# dataに含まれるデータ数の確認 p.15
len(data)

【実行結果】
50 個（50 行）のデータです。

③ 要約統計量の表示
describe メソッドで要約統計量を表示できます。
１列しかないので、横並びに表示すると見やすいです。
T 属性で縦横変換（行列の転置に相当）ができます。 

### dataの要約 要約統計量の表示 p.16
data.describe().T

【実行結果】
左からデータの個数、平均、標準偏差、最小値、第１四分位数、中央値、第三四分位数、最大値です。
標準偏差は不偏分散の標準偏差です。

平均 3.56 を覚えておきます。後で分散と比べます。

④ 度数分布の表示
value_counts メソッドで要素の数（ユニークな種子数）と出現頻度を取得します。

### 度数分布を表示 p.16
data.value_counts(sort=False).to_frame().T

【実行結果】
種子数は０以上の整数値です。
種子数は 0 から 7 の範囲にあり、度数の最大値は 12（種子数 3 のとき）です。

⑤ ヒストグラムの描画
pandas の plot メソッドでグラフを描きましょう。
データフレーム data に対する plot.hist メソッドでヒストグラムを描画できます。

### ヒストグラムの描画 p.17 図2.2
data.plot.hist(
    bins=np.arange(-0.5, 8), xticks=range(9), edgecolor='white', alpha=0.7,
    legend=False, title='Histogram of data', xlabel='種子数'
);

【実行結果】
峰（山）が一つの分布です。

⑥ 標本分散・標本標準偏差の算出
分母を $${n-1}$$（分散の場合は不偏分散）にして計算する場合、引数 ddof を 1 にします。 
今回使用するメソッドの場合、ddof のデフォルト値が 1 ですので、コードでは設定を省略しています。

標本分散

### 標本分散の表示 データのばらつき p.17
data['種子数'].var()  # ddof=1

【実行結果】
標本分散 2.99 は標本平均 3.56 と近い値です（要記憶）。

標本標準偏差

### 標本標準偏差の表示 p.17
data['種子数'].std()  # ddof=1

【実行結果】

標本標準偏差（標本分散の正の平方根）

### 標本標準偏差の表示 ※標本分散の平方根 p.17
np.sqrt(data['種子数'].var())  # ddof=1

【実行結果】

◆ ◆ ◆

■ データの特徴まとめ
テキスト p.18 冒頭のまとめと p.22 の確率分布選定理由をお借りして、データの特徴を整理します。

データのばらつきは「確率分布」で表現します。
上述のデータの特徴を表現できる確率分布は「ポアソン分布」です。

統計モデルをデータに当てはめ（モデルの理解）

今回はポアソン分布を例題データ＝観測データに当てはめます。

■ ポアソン分布の概要
ポアソン分布は確率変数がとびとびの値（離散値）をとる離散型確率分布です。
確率変数を $${y}$$、平均パラメータを $${\lambda}$$ とすると、ポアソン分布の確率質量関数は次の式で表されます。

$${!}$$ は階乗の記号です。
たとえば $${y=3}$$ のとき、$${3! = 3 \times 2 \times 1 = 6}$$ です。

ポアソン分布のプロフィールです。

Python でポアソン分布の確率質量関数を表現します。
① 確率質量関数の数式を関数化

# ポアソン分布の確率質量関数を算出する関数の定義
def poisson_pmf(y, lam):
    return lam**y * np.exp(-lam) / np.array([math.factorial(i) for i in y])

確率変数の範囲 $${y = 0, 1, \cdots ,9}$$、平均パラメータ $${\lambda = 3.5}$$ で確率質量関数を計算します。

# データyと平均パラメータλの設定
y = np.arange(10)
lam = 3.5

# 確率の算出
poisson_pmf(y=y, lam=lam)

【実行結果】

② scipy.stats の ポアソン分布クラス利用

# scipy.statsのポアソン分布の確率質量関数

# 確率の算出
stats.poisson.pmf(k=y, mu=lam)

【実行結果】

ポアソン分布の確率質量関数を可視化します。
テキスト p.23 の図 2.6 に相当します。
確率変数の範囲 $${y = 0, 1, \cdots ,21}$$、平均パラメータ $${\lambda = \{3.5, 7.7, 15.1\}}$$ のケースで描画します。

### さまざまな平均(λ)のポアソン分布 p.23 図2.6

## 設定と準備
lams = [3.5, 7.7, 15.1]                        # 平均λ
markers = ['o', 'D', '^']                      # マーカーの形状
colors = ['tab:blue', 'tab:red', 'tab:green']  # グラフの色
y = np.arange(21)                              # x軸の値（確率変数yの値）

## 描画
# ３つの平均ごとに確率算出と折れ線グラフ描画を繰り返し処理
for lam, marker, color in zip(lams, markers, colors):
    # 平均λを用いて確率を算出
    prob = stats.poisson.pmf(k=y, mu=lam)
    # 折れ線グラフの描画
    plt.plot(y, prob, c=color, mfc='white', marker=marker, ms=8, ls='--',
             label=f'{lam}')
# 修飾：x軸ラベル、y軸ラベル、x軸目盛り、凡例
plt.xlabel('y', fontsize=12)
plt.ylabel('prob', fontsize=12)
plt.xticks(ticks=range(0, 21, 5))
plt.legend(title='lambda');

【実行結果】
分布の形状は峰が一つです。
$${\lambda}$$ の値が大きくなるにつれて、峰が上側（右側）に移動し、尖り具合が潰れてバラツキが大きくなります。

◆ ◆ ◆

■ 観測データとポアソン分布
観測データの特徴はポアソン分布とよく似ています。

ポアソン分布を観測データに当てはめる際の「前提条件」をテキストから引用いたします。

統計モデルをデータに当てはめ（当てはめと評価）

まずテキスト 2.2 節の「データとポアソン分布の対応」を実践します。
次にテキスト 2.4 節の「最尤推定」を実践します。

■ $${\lambda =}$$ データの平均値のポアソン分布を当てはめる
観測データの平均値 $${3.56}$$ をポアソン分布の平均パラメータ $${\lambda}$$ に設定して、データとポアソン分布の対応を確認します。

確率変数である種子数 $${y}$$ が $${0, 1, 2, \cdots, 9}$$ のときのポアソン分布の確率質量関数を計算します。
テキスト p.19 の最初の R コードに対応します。
scipy.stats の poisson クラスを利用します。

### 平均3.56のポアソン分布の確率分布の計算

# 確率変数yの値の設定
y = np.arange(10)

# yについて平均3.56のポアソン分布の確率質量関数を算出 scipy.stats利用
prob = stats.poisson.pmf(k=y, mu=data.mean())
prob

【実行結果】

p.20 図2.3 に相当する確率質量関数の表を作成します。
（上の prob を pandas データフレームに設定するだけです）

# データフレーム化 p.20 図2.3
prob_df = pd.DataFrame({'y': y, 'prob': prob})
prob_df

【実行結果】
種子数 3 のときに確率の最大値になっています。
平均値がピークになっている感じです。

この確率値を可視化しましょう。
$${\lambda=3.56}$$ のポアソン分布の確率質量関数を描きます。
テキスト p.20 図 2.4 に相当します。
pandas の plot メソッドを使います。

### yと確率の描画 p.19, 20 図2.4
prob_df.plot(x='y', y='prob', marker='o', mfc='white', ls='--',
             ylabel='prob', legend=False);

【実行結果】
先の表の通り、種子数 3 のときが確率の峰（ピーク）になっていて、少々下側（左側）に偏った分布になっています。

いよいよ最初のクライマックスです！
「観測データの分布を示すヒストグラム」と「ポアソン分布による種子数の予測値」（確率質量関数 $${\times}$$ 全個体数 $${50}$$）を重ね描きしましょう。
テキスト p.21 図 2.5 に相当します。
pandas の plot メソッドを使います。

### yと確率の描画 p.21 図2.5

# ヒストグラムの描画
ax = data.plot.hist(bins=np.arange(-0.5, 8), xticks=range(9),
                    color='tab:blue', edgecolor='white', alpha=0.7,
                    legend=False, title='Histogram of data')
# 個体数の予測値の折れ線グラフの描画
prob_df['pred'] = prob_df['prob'] * len(data)
prob_df.plot(x='y', y='pred', marker='o', color='tab:red', mfc='white',
             lw=3, ls='--', legend=False, ax=ax);

【実行結果】

【考察】
観測データのヒストグラムの分布（青い棒）に対して $${\lambda=3.56}$$ のポアソン分布による予測値（赤い点線）が「よく当てはまっている」感じがいたします！

これが「ポアソン分布を観測データに当てはめる」ということでしょうか！？
ポアソン分布のパラメータ $${\lambda}$$ を観測データから推定することが「ポアソン分布をデータに当てはめる」ことなのでしょうか！？（そうです🌛）

テキストの言葉をお借りしますと…

定量的に示す方法、つまり、確率分布のパラメータを定量的に決める方法の一つが「最尤推定」です。

（注）
テキスト p.21 の「観測データを説明できる良い統計モデルとは何か」のアンサーは「予測の良さ」（新しく得られた未知データへの当てはまりの良さ）と思われます。

◆ ◆ ◆

■ 最尤推定で平均パラメータ $${\lambda}$$ を推定する
観測データに基づいて確率分布のパラメータを推定します。
テキストは「最尤推定」でパラメータ推定を行います。

最尤推定は、観測データから「尤度」と呼ばれる統計量を最大にするパラメータの値を推定する方法です。

① 尤度関数と対数尤度関数
尤度は、パラメータ $${\theta}$$ をある値に定めたとき、全ての観測データ $${y_i}$$ についての確率 $${p(y_i \mid \theta)}$$ の積です。

📊 尤度関数の公式

$${\bm Y}$$ は観測値（ベクトル）、$${N}$$ は標本サイズ、$${\prod}$$ は総乗記号です。

尤度関数を対数変換した「対数尤度関数」を一般に用いるようです。
確率値の掛け算は非常に小さな値になるため避けたい、とか、「積」を「和」にすることで微分計算がしやすくなることなどが理由だそうです。

📊 対数尤度関数の公式

② ポアソン分布の尤度関数・対数尤度関数
ポアソン分布の確率質量関数：

を用いて、尤度関数は次のようになります。

尤度関数を対数変換した対数尤度関数は次のようになります。

観測データに関する対数尤度関数を Python で実装します。
上述の数式を関数化して、$${\lambda=3.6}$$ のときの対数尤度関数を計算します。

# ポアソン分布の対数尤度関数 ※p.26の数式利用
def logL(λ, y):
    return sum(yi * np.log(λ) - λ - sum([np.log(k) for k in range(1, yi+1)])
               for yi in y)

logL(λ=3.6, y=data['種子数'])

【実行結果】
対数尤度関数 $${\log L(\lambda=3.6)}$$ は $${-97.3}$$ です。

scipy.stats で対数尤度関数を実装します。
観測データ１点ごとに対数確率質量関数を計算して合計します。

# 次の図で用いるポアソン分布の対数尤度関数 ※scipy.stats利用
np.sum(stats.poisson.logpmf(k=data['種子数'], mu=3.6))

【実行結果】
対数尤度関数 $${\log L(\lambda=3.6)}$$ は $${-97.3}$$ です。

◆ ◆ ◆

③ ポアソン分布を観測データに当てはめる実験
パラメータ $${\lambda}$$ の値を変化させることで、ポアソン分布の確率質量関数を観測データの分布に近づける実験をテキストが行っています。
Python でこの実験を実装します。
p.26 の図 2.7 に相当します。

### 平均λを変化させてポアソン分布と観測データの当てはまりの良さを調べる p.27 図2.7

## 設定
lams = np.arange(2.0, 5.3, 0.4)  # 平均λ
y = range(10)                    # x軸の値（確率変数yの値）

## 描画
# 描画領域の設定
fig, axes = plt.subplots(3, 3, figsize=(15, 12), tight_layout=True,
                         sharex=True, sharey=True)

# 平均λの値ごとにポアソン分布の確率算出・対数尤度算出・チャート描画を繰り返し処理
for lam, ax in zip(lams, axes.flat):

    ## 描画・表示用データの算出
    # 平均λを用いてポアソン分布の確率質量関数を算出しデータフレーム化
    prob = stats.poisson.pmf(k=y, mu=lam)
    prob_df2 = pd.DataFrame({'y': y, 'prob': prob})
    # 種子数データを用いて対数尤度の算出
    log_L = np.sum(stats.poisson.logpmf(k=data['種子数'], mu=lam))

    ## 描画
    # 種子数のヒストグラムの描画
    data.plot.hist(bins=np.arange(-0.5, 8), xticks=range(9), ylim=(0, 15),
                   color='tab:blue', edgecolor='white', alpha=0.5,
                   legend=False, ax=ax)
    # ポアソン分布の確率質量関数の描画
    ax2 = ax.twinx()  # y軸目盛りを変更（右側に表示）
    prob_df2.plot(x='y', y='prob', color='tab:red', ylim=(0, 0.3), 
                  marker='o', ms=10, mfc='white', lw=3, ls='--', legend=False,
                  ax=ax2)
    ax2.set_yticks([])
    # 平均λと対数尤度の表示
    ax.text(x=6, y=12, s=f'lambda={lam:.1f}\nlog L={log_L:.1f}', fontsize=14)

【実行結果】
パラメータ $${\lambda}$$ の値を変えてポアソン分布を観測データに当てはめ（フィッティング）するイメージが掴めるような気がします。

【考察】
対数尤度 $${\log L(\lambda \mid \bm Y)}$$ が最大になる $${\lambda}$$ は $${3.2}$$ と $${3.6}$$ の間にありそうです。
$${\lambda=3.6}$$ のとき、観測データのヒストグラムの形状とポアソン分布の確率質量関数の形状が最も近くなっている＝良く当てはまっているように見えます。

◆ ◆ ◆

④ 最尤推定の実行
先ほどのポアソン分布の当てはめ実験では、$${\lambda}$$ の値を $${0.4}$$ 間隔にしました。
さらに間隔を小さくして、対数尤度が最大となるパラメータ $${\lambda}$$ の値を探索しましょう。
テキスト p.28 図 2.8 に相当します。

### 種子数データのもとでのλと対数尤度の関係 p.28 図2.8

## 描画データの作成
# 対数尤度算出関数の定義 scipy.stats利用
logL = lambda lams, y: np.array(
    [np.sum(stats.poisson.logpmf(k=y, mu=lam)) for lam in lams])
# x軸の平均パラメータλの値
lams = np.linspace(2, 5, 1001)
# 平均パラメータλに対応する対数尤度の算出
log_Ls = logL(lams, y=data['種子数'])
# 対数尤度の最大値のインデックスの取得
max_idx = log_Ls.argmax()

## 描画
# λと対数尤度の曲線の描画
plt.plot(lams, log_Ls)
# 対数尤度が最大となるλの垂直線の描画
plt.axvline(lams[max_idx], color='black', ls='--', lw=0.9)
# 曲線と垂直線の接点の描画
plt.plot(lams[max_idx], log_Ls[max_idx], color='tab:red', marker='x', ms=14)
# 対数尤度が最大となるλの値の表示
plt.text(x=3.7, y=-122, s=f'$\hat{{\lambda}}$={lams[max_idx]:.2f}', fontsize=14)
# 修飾：x軸ラベル、y軸ラベル、タイトル
plt.xlabel('$\lambda$', fontsize=14)
plt.ylabel('対数尤度', fontsize=14)
plt.title(f'最大対数尤度 {log_Ls[max_idx]:.3f}');

【実行結果】
横軸が $${\lambda}$$、縦軸が $${\lambda}$$ に対応する対数尤度  $${\log L(\lambda \mid \bm Y)}$$ です。

【最尤推定の結果】
対数尤度が最大となる $${\lambda}$$ は $${3.56}$$ であり、最大対数尤度は $${-97.24}$$ になりました。
パラメータ $${\lambda}$$ の最尤推定値は $${\hat{\lambda} = 3.56}$$ です。
最尤推定値は観測データの標本平均と同じになりました。

◆ ◆ ◆

■ 解析的に最尤推定量 $${\hat{\lambda}}$$ を計算する
テキスト p.27 では解析的に最尤推定量を算出しています。
対数尤度関数をパラメータ $${\lambda}$$ で偏微分して、偏微分結果をゼロとおいて求めます。

【偏微分】

【偏微分＝０を求解】

【結論】
ポアソン分布のパラメータ $${\lambda}$$ の最尤推定量 $${\hat{\lambda}}$$ は標本平均です。
観測データの場合、最尤推定値は $${\hat{\lambda} = 3.56}$$ です。

◆ ◆ ◆

■ 閑話休題：最適化ライブラリで最尤推定
scipy の 最小化問題ソルバー minimize を用いて最尤推定をやってみます。
「対数尤度最大化」を「負の対数尤度の最小化」に変えるところがミソです。

# 最適化ライブラリで最尤推定

# 追加インポート
from scipy.optimize import minimize

# 対数尤度関数（目的関数）※負の対数尤度の最小化問題を解く
fun = lambda lam, x: -np.sum(stats.poisson.logpmf(k=x, mu=lam))

# λの最尤推定値を算出
res_min = minimize(fun=fun, x0=[0.01], args=(data['種子数']))

# 結果の表示
print(f'λの最尤推定値 = {res_min.x[0]:7.3f}')
print(f'最大対数尤度  = {-res_min.fun:7.3f}')

【実行結果】
とてもお手軽に最尤推定できました！

◆ ◆ ◆

■ 最尤推定値のバラツキの評価
最尤推定値のバラツキ＝標準誤差 $${\text{SE}}$$ を見積もります。
テキスト p.29～の乱数生成による見積もり方法を実践します。
出来上がりの図はテキスト p.30 図2.9 に相当します。

ここでは「パラメータ $${\lambda}$$ の真値を知っている」ことを前提とします。
$${\lambda=3.5}$$ のポアソン分布から 50 個乱数を生成して最尤推定値 $${\hat{\lambda}}$$ ＝標本平均を計算します。
この最尤推定値を 3000 個分採取します。
3000 個のデータの標本標準偏差が 最尤推定値 $${\hat{\lambda}}$$ の標準誤差です。

### ポアソン分布の平均λの最尤推定値λ_hatのばらつき p.30 図2.9

# 設定と準備
n = 50               # サンプルサイズ
lam_true = 3.5       # ポアソン分布の平均λ
iter = 3000          # シミュレーション回数
rng = np.random.default_rng(seed=123)  # 乱数生成器
samples = []         # λ_hatを格納するリストの初期化

# λの最尤推定を3000回繰り返す
for _ in range(iter):  # ポアソン分布のλの最尤推定値が標本平均であることを利用
    lam_hat = stats.poisson.rvs(mu=lam_true, size=n, random_state=rng).mean()
    samples.append(lam_hat)

# ヒストグラムの描画
plt.hist(samples, bins=20, edgecolor='white', alpha=0.7)
plt.xlabel('試行ごとに推定された $\hat{\lambda}$', fontsize=12);

【実行結果】
標本平均 $${3.56}$$ を中心にして、$${2.6}$$ あたりから $${4.3}$$ くらいまでばらついています。

最尤推定値の標準誤差を計算します。
3000 個のサンプルの標準偏差です。

# シミュレーションに基づく最尤推定値の標準誤差
print(f'最尤推定値の標準誤差: {np.std(samples, ddof=1):.4f}')

【実行結果】
標準誤差は約 $${0.27}$$ です。

標準誤差 $${\text{SE}}$$ の理論値は：

となります。

テキスト p.30 に「推定値のばらつきを見積もるときには、$${\hat{\lambda}}$$ を利用」するとの記述があります。

$${\hat{\lambda} = 3.56}$$ を用いて標準誤差を見積もってみます。

### ポアソン分布の平均λの最尤推定値λ_hatのばらつき 最尤推定値で計算

# 設定と準備
n = 50               # サンプルサイズ
lam_llf = 3.56       # ポアソン分布の平均λの最尤推定値
iter = 3000          # シミュレーション回数
rng = np.random.default_rng(seed=123)  # 乱数生成器
samples2 = []         # λ_hatを格納するリストの初期化

# λの最尤推定を3000回繰り返す
for _ in range(iter):  # ポアソン分布のλの最尤推定値が標本平均であることを利用
    lam_hat = stats.poisson.rvs(mu=lam_llf, size=n, random_state=rng).mean()
    samples2.append(lam_hat)

# ヒストグラムの描画
plt.hist(samples2, bins=20, edgecolor='white', alpha=0.7)
plt.xlabel('試行ごとに推定された $\hat{\lambda}$', fontsize=12);

【実行結果】

# シミュレーションに基づく最尤推定値の標準誤差
print(f'最尤推定値の標準誤差: {np.std(samples2, ddof=1):.4f}')

【実行結果】
真値 $${3.5}$$ を用いたときよりも標準誤差は（若干）大きくなりました。

◆ ◆ ◆

■ 統計モデルの予測の良さの評価
４章で学ぶ「AIC」を先取りして計算します。
AIC は複数モデルから「もっとも予測の良いモデル」を選択する際に使用する指標であり、AICが小さいモデルほど「予測の良いモデル」とされます。

📊 AIC の公式

$${\log L^*}$$ は最大対数尤度、$${k}$$ はパラメータの数です。
今回のポアソン分布モデルのパラメータ数は１です。

# 統計モデルのモデルの予測の悪さ AIC の算出

# 最大対数尤度の取得
llf_max = -res_min.fun

# AICの算出 ※AIC = -2 * (最大対数尤度 - パラメータ数) 
aic = -2 * (llf_max - 1)
print(f'AIC = {aic:.3f}')

【実行結果】

予測

テキストでは p.21 図 2.5 の赤い点線「種子数の予測値」が予測にあたります。
最尤推定値 $${\hat{\lambda} = 3.56}$$ のポアソン分布の確率質量関数の値（図 2.3 の表の値）に全個体数 50 を乗じて計算しています。

アディショナル：GLMの練習

今回のポアソン分布の当てはめを一般化線形モデル GLM で表現してみます。
３章に入る前の GLM の練習みたいなものです。

■ 確率分布、リンク関数、線形予測子
GLMの３要素である「確率分布」、「リンク関数」、「線形予測子」を定義します。

🔷 確率分布と確率質量関数
種子数 $${y_i}$$ はパラメータ $${\lambda}$$ のポアソン分布に従います。

🔷 リンク関数と線形予測子
リンク関数は「対数」、線形予測子は「切片 $${\beta_0}$$」です。

◆ ◆ ◆

■ GLM の実行
statsmodels の glm を利用して上記の統計モデルを実装します。

# GLMで最尤推定

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

# 確率分布＝ポアソン分布、リンク関数＝Logの設定（リンク関数の設定は省略可能）
family=sm.families.Poisson(link=sm.families.links.Log())

# 線形予測子をformulaで指定
formula = '種子数 ~ 1'  # 定数項のみのモデル

# statsmodelsのglm関数を利用してポアソン回帰を実行
result = smf.glm(formula=formula, data=data, family=family).fit()
result.summary()

【実行結果】

【考察】
対数尤度は $${-97.244}$$ です。
いままで検討したポアソン分布単独のモデルの最大対数尤度と同じ値になっています。

◆ ◆ ◆

■ パラメータ $${\lambda}$$ の推定
GLMの結果では、ポアソン分布のパラメータ $${\lambda}$$ は明示されません。
線形予測子の $${\beta_0}$$（切片）が推定されていますので、$${\lambda = \exp(\beta_0)}$$ を用いてパラメータ $${\lambda}$$ を推定します。

# ポアソン分布の平均パラメータλの算出
lam = np.exp(result.params['Intercept'])
print(f'λ = {lam:.4f}')

【実行結果】
いままで検討したポアソン分布単独のモデルの最尤推定値と同じ値になっています。

■ $${\hat{\lambda}}$$ の標準誤差の見積もり
最尤推定値 $${\hat{\lambda}}$$ の標準誤差の見積もりは少々難しいです。
ChatGPT に相談したところ、「delta法」で標準誤差の近似値を求められるとのこと。
早速やってみます。

# λの標準誤差を推定する

# ポアソン回帰の結果から切片の標準誤差を取得
se_beta0 = result.bse['Intercept']

# delta法を用いて平均パラメータλの標準誤差を算出
se_lam = lam * se_beta0

# 結果の表示
print(f'標準誤差 = {se_lam:.4f}')

【実行結果】
ポアソン乱数生成で推定した標準誤差 $${2.668}$$ と一致しています。

■ AIC の算出
GLM の実行結果から AIC を取り出します。

# AICの取得
print(f'AIC = {result.aic:.3f}')

【実行結果】
いままで検討したポアソン分布単独のモデルの AIC と同じ値になっています。

◆ ◆ ◆

■ 予測
GLM の結果を用いて予測をします。
matplotlib の plot 関数を利用して可視化します。

# 予測

# 予測値の算出
y_pred = result.predict()

# 観測データの散布図の描画
plt.plot(data, ls='none', marker='o', alpha=0.7, label='観測値')
# 予測値の直線の描画
plt.plot(y_pred, ls='--', color='tab:red', label='予測値')
# 修飾
plt.xlabel('データ番号', fontsize=12)
plt.ylabel('種子数', fontsize=12)
plt.legend();

【実行結果】
青い点の観測値はばらついています。
赤い点線の予測値は $${3.56}$$ で一定になっています。

【考察】
50 個の観測データに対応する予測値はすべて $${3.56}$$（$${\lambda}$$ の最尤推定値）になっています。
切片モデルなので、こうなってしまうのでしょう。

その他のテキストの見どころ

■ 統計モデルの要点
データ生成（乱数生成）⇒パラメータ推定⇒モデルによる予測の一連の流れ・考え方は、テキスト p.32 図 2.10、p.33 図 2.11 をぜひご確認ください。

■ 確率分布の選び方 
確率分布の選び方のガイドラインをテキスト p.34 より引用いたします。

■ 観測データの分布が単純な確率分布に見えない場合 p.35
「データ化されていない個体差・観測しなかった個体差」を組み込んだ統計モデリングによる対応をテキストは提案しています。
おそらく 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の教科書です。
よかったらぜひ、お試しくださいませ。

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