StanとRでベイズ統計モデリングをPyMC Ver.5で写経～第10章「10.1.4 多項ロジスティック回帰」

2024年3月19日 11:44

第10章「収束しない場合の対処法」

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

この記事は、テキスト第10章「収束しない場合の対処法」・10.1節「パラメータの識別可能性」の 10.1.4項「多項ロジスティック回帰」の PyMC5写経 を取り扱います。
実は初回モデリング時には収束しませんでした（笑）
モデルのパラメータ値を試行錯誤で求める、という対処法で収束に漕ぎ着けました。

テキストは第10章で収束に向けた工夫を取り扱っています。
座学や数式モデルのみの掲載項があれば、Stanコードの掲載項もあります。
PyMC化は主にStanコードが明示されているモデル式を対象にして実施します。

はじめに

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を動かすまでの準備」章をご覧ください。

10.1.4 多項ロジスティック回帰

インポート

### インポート

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

# PyMC
import pymc as pm
import pytensor.tensor as pt
import arviz as az

# 評価指標
from sklearn.metrics import classification_report

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

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

データの読み込み・確認

サンプルコードのデータを読み込みます。

### データの読み込み ◆データファイル10.1 data-category.txt
# 購入者属性 Age:年齢, Sex:性別(0:男性,1:女性), Income:年収[万円]
# 目的変数   Y:購入商品のカテゴリーのインデックス(1～Kのいずれかの整数値)

data = pd.read_csv('./data/data-category.txt')
print('data.shape: ', data.shape)
display(data.head())

【実行結果】

集約・可視化を行ってみます。
まずは集約・可視化しやすいようにデータを加工します。

### 集約・可視化に用いるデータの作成
# dataのコピーを取得
sum_data = data.copy()
# 所得層項目の追加（５段階）
sum_data['Income_label'] = pd.cut(data['Income'], 5,
                    labels=['x-small', 'small', 'medium', 'large', 'x-large'])
# 年代項目の追加（~29、~39, ~49, 50~）
sum_data['Ages'] = (sum_data['Age'].apply(lambda x:
                    20 if x<30 else 30 if x<40 else 40 if x<50 else 50))

【実行結果】なし

年代別・商品カテゴリ別の購入数を集計します。

### 年代別・商品カテゴリ別の購入数
sum_data.pivot_table(index='Ages', columns='Y', values='Age', aggfunc='count')

【実行結果】

この集計データを可視化してみます。

### 年代別・商品カテゴリ別の購入数の描画
sns.countplot(data=sum_data, x='Y', hue='Ages', palette='tab10', alpha=0.7)
plt.legend(bbox_to_anchor=(1.18, 1), title='年代')
plt.grid(lw=0.5);

【実行結果】
年代によって購入数の多い／少ないが見えてきましたか？

次は性別・商品カテゴリ別の購入数を集計します。

### 性別・商品カテゴリ別の購入数
sum_data.pivot_table(index='Sex', columns='Y', values='Age', aggfunc='count')

【実行結果】
性別は、0: 男性、1: 女性です。

この集計データを可視化してみます。

### 性別・商品カテゴリ別の購入数の描画
sns.countplot(data=sum_data, x='Y', hue='Sex', palette='tab10', alpha=0.7)
plt.legend(bbox_to_anchor=(1.18, 1), title='性別', labels=['男性', '女性'])
plt.grid(lw=0.5);

【実行結果】
男性と女性の購入数の違いが見えましたね！

最後に所得層別・商品カテゴリ別の購入数を集計します。

### 所得層別・商品カテゴリ別の購入数
sum_data.pivot_table(index='Income_label', columns='Y', values='Age',
                     aggfunc='count')

【実行結果】

この集計データを可視化してみます。

### 所得層別・商品カテゴリ別の購入数の描画
sns.countplot(data=sum_data, x='Y', hue='Income_label', palette='tab10',
              alpha=0.7)
plt.legend(bbox_to_anchor=(1.25, 1), title='所得層')
plt.grid(lw=0.5);

【実行結果】
購買層は見えてきましたでしょうか？

PyMCのモデル定義

PyMCでモデル式10-2を実装します。
まずは説明変数$${X}$$を整えます。

### 説明変数 x の作成
x = data.iloc[:, :3]
x['const'] = 1
x = x[['const', 'Age', 'Sex', 'Income']]
display(x)

【実行結果】

モデルの定義です。
テキストに掲載の識別不可能を回避する対策「ベクトル$${\overrightarrow{b_1}}$$の１番目の要素を０に固定すればよい」を b0 で表現しています。
また収束のための試行錯誤は、b1 の一様分布の範囲をアレコレ試して、$${[-1, 1]}$$で収束できました。

### モデルの定義 ◆モデル式10-2 model10-2.stan

with pm.Model() as model:
    
    ### データ関連定義
    ## coordの定義
    # データのインデックス
    model.add_coord('data', values=data.index, mutable=True)
    # Yのカテゴリ １～６
    model.add_coord('category', values=sorted(data['Y'].unique()), mutable=True)
    # ４つの回帰係数
    model.add_coord('coef', values=x.columns, mutable=True)
    ## dataの定義
    # 目的変数「Y - 1」 0始まり化
    Y = pm.ConstantData('Y', value=data['Y'].values - 1, dims='data')
    # 説明変数 X
    X = pm.ConstantData('X', value=x.values, dims=('data', 'coef'))

    ### 事前分布
    ## 回帰係数 b
    # カテゴリ１を０に固定（4行1列）
    b0 = pt.zeros((4, 1))
    # カテゴリ２～６（4行5列）
    b1 = pm.Uniform('b1', lower=-1, upper=1, shape=(4, 5))
    # カテゴリ１とカテゴリ２～６を列方向に結合
    b = pm.Deterministic('b', pt.concatenate([b0, b1], axis=1),
                         dims=('coef', 'category'))
    ## 線形予測子 μ
    mu = pm.Deterministic('mu', pt.dot(X, b), dims=('data', 'category'))
    ## 確率 θ （softmax関数適用）
    theta = pm.Deterministic('theta', pm.math.softmax(mu, axis=1),
                             dims=('data', 'category'))

    ### 尤度関数
    obs = pm.Categorical('obs', p=theta, observed=Y, dims='data')

モデルの定義内容を見ます。

### モデルの表示
model

【実行結果】

### モデルの可視化
pm.model_to_graphviz(model)

【実行結果】

MCMCの実行と収束確認

MCMCを実行します。

### 事後分布からのサンプリング 50秒
with model:
    idata  = pm.sample(draws=1000, tune=1000, chains=4, target_accept=0.8,
                       nuts_sampler='numpyro', random_seed=1234)

【実行結果】省略

Pythonで事後分布からのサンプリングデータの確認を行います。
Rhatの確認から。
テキストの収束条件は「chainを3以上にして$${\hat{R}<1.1}$$のとき」です。

### r_hat>1.1の確認
# 設定
idata_in = idata         # idata名
threshold = 1.01         # しきい値

# しきい値を超えるR_hatの個数を表示
print((az.rhat(idata_in) > threshold).sum())

【実行結果】
収束条件を満たしています。

事後統計量を表示します。

### 推論データの要約統計情報の表示
var_names = ['b', 'mu', 'theta']
pm.summary(idata, hdi_prob=0.95, var_names=var_names, round_to=3)

【実行結果】

トレースプロットを描画します。
パラメータの一部を取り上げています。

### トレースプロットの表示
pm.plot_trace(idata, var_names=var_names, compact=False)
plt.tight_layout();

【実行結果】

Y の事後予測サンプリングを実行します。

### Yの事後予測サンプリングの実行
with model:
    pp = pm.sample_posterior_predictive(idata, random_seed=1234)

【実行結果】

事後予測プロットでモデルの予測の良し悪しを確認しましょう。

### 事後予測プロットの描画
pm.plot_ppc(pp, num_pp_samples=100, random_seed=1234);

【実行結果】
観測値と事後予測の平均値は一致している感じですが、個々のサンプリングデータ（青線）は振れ幅が大きい感じもします。

推定結果の解釈

多値分類予測の評価指標を見てみます。
確率$${\theta}$$の平均値が最も大きいカテゴリを予測値に見立てました。

### 予測精度を見てみる

# Yの予測値 theta(確率値)の平均値が最も大きいカテゴリを予測値にする
y_pred = np.array([(idata.posterior.theta.stack(sample=('chain', 'draw'))[i]
                    .mean(axis=1).argmax().values + 1)
                    for i in range(300)])

# Yの観測値（正解値）
y_true = data['Y'].values

# 分類レポートの表示
print(classification_report(y_true, y_pred))

【実行結果】
正解率（accuracy）は$${0.31}$$。精度は低そうです。
上段を見ますと、カテゴリ３と５は全く当てられていません。

10.1.4 項は以上です。

シリーズの記事

次の記事

前の記事

ブログの紹介

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

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

いいなと思ったら応援しよう！

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

#うちの積読を紹介する

2,210件