見出し画像

第14章「歯医者さんに通うのはお好きですか？」のベイズモデリングをPyMC Ver.5 で

ネイピア DS

この記事は、テキスト「たのしいベイズモデリング２」の第14章「歯医者さんに通うのはお好きですか？」のベイズモデルを用いて、PyMC Ver.5で「実験的」に実装する様子を描いた統計ドキュメンタリーです。

章タイトルのインパクトが強いですね！

この章では、歯科医と患者との間に「価値共創」を生み出す患者の行動をベイズ推論します。
初学者にやさしい線形回帰モデルでマーケティング論に軽くタッチできる章です。

歯医者のイラスト「治療中」：「いらすとや」さんより

今回の PyMC化 は、テキストの結果と概ね一致しました！
ではでは、PyMCのベイズモデリングの世界を楽しんでまいりましょう！

テキストの紹介、引用表記、シリーズまえがき、PyMC等のバージョン情報は、このリンクの記事をご参照ください。

テキストで使用するデータは、R・Stan等のサンプルスクリプトとともに、出版社サイトからダウンロードして取得できます。


サマリー

テキストの概要

執筆者　　　：　五島光 先生
モデル難易度：　★・・・・　（やさしい）

自己評価

評点

$$
\begin{array}{c:c:c}
実装精度 & ★★★★★ & GoooD! \\
結果再現度 & ★★★★★ & 最高✨ \\
楽しさ & ★★★★★ & 楽しい！ \\
\end{array}
$$

花型の評価印のイラスト：「いらすとや」さんより

評価ポイント

  • テキストの結果とほぼ一致しました。
    調査目的指向のシンプルなモデリングはいいですね～

工夫・喜び・反省

  • テキストのおわりの一文「本章で扱ったように Stan とそのラッパーを使えば、マーケティングサイエンスの専門家ではない統計学エンドユーザーにもベイズ推定が手に届くようになってきた」（太字は私が付けました）に注目です！
    PyMC にも心強い親友 Bambi が居ます！
    Bambi を使えば専門家ではない統計学エンドユーザーにベイズ推定が手に届くと思います！

子鹿・バンビのイラスト：「いらすとや」さんより

モデルの概要

テキストの調査・実験の概要

■ 価値共創
この章の分析の目的をテキストから引用いたします。

私たちの生活を向上させるサービス（歯科医院）において、よりよい価値共創が起きるためにサービスの受け手がどのような関わり方をすればいいのかを明らかにする。

テキストより引用

サービスを提供する「歯科医院」とサービスを受け取る「患者」の相互関わり合い（インタラクション）によって成果が生み出されること～価値共創～。
価値共創が起きるために患者サイドがどのような関わりをするのがいいのでしょう？

◆ ◆ ◆ ◆ ◆

■ 価値共創実践スタイル（CVCPS）
テキストは、サービスの価値共創における受け手の関わり方として「価値共創実践スタイル（CVCPS）」をベースにして、議論を展開します。
歯科医院における５つのCVCPSを提案して分析を進めます。

【歯科医院におけるCVSPS】
v1：助言・指示の遵守：患者が歯科医院の指示を守ること
v2：積極的コミットメント：患者が積極的に指示や説明を求めること
v3：治療選択：患者自ら考えて判断すること
v4：私的治療等会話：患者がプライベートで口・歯を話題にすること
v5：健康情報意識：患者が口・歯の健康情報に注目すること

テキストの記述を一部改変して引用
歯医者さんのイラスト（女性）：「いらすとや」さんより

■ データの概要
１年以内に歯科医院に通院した人を対象にしたオンラインリサーチの方法で、上記のCVCPS５項目と、歯科医院に通うことで生活の質が高まったか（目的変数：QOL）を１～７段階（１が最も低く７が最も高い）で質問し、500人から回答を集めたそうです。
質問文の詳細はテキストでご確認ください。

ここまでの概要を図にしました。

テキストのモデリング

■ 目的変数と関心のあるパラメータ
QOLの回答値$${QOL}$$を目的変数、CVCPS５項目の回答値$${v_1}$$～$${v_5}$$を説明変数にする線形回帰モデルを構築します。
注目するパラメータはCVCPS５変数の（偏）回帰係数です。

■ モデル数式
テキストは brms を利用しています。
brms は「Bayesian Regression Models using 'Stan'」（Stanを利用したベイジアン回帰モデル）の略で、R の Formula 構文でモデリングできます。
Formula 構文による今回のモデルは次のとおりです。

$$
QOL \sim v_1 + v_2 + v_3 + v_4 + v_5 \\
$$

Stanスクリプトの記述を一部改変して引用

確率変数を用いる記述に変換するとおそらくこんな感じです。

$$
\begin{align*}
QOL &\sim \text{Normal}\ (\mu,\ \sigma^2) \\
\mu &= \beta_0 + \beta_1v_1 + \beta_2v_2 + \beta_3v_4 + \beta_4v_4 + \beta_tv_5 \\
\end{align*}
$$

■分析・分析結果
分析の仕方や分析数値はテキストの記述が正確だと思いますので、テキストの読み込みをおすすめします。
PyMCの自己流モデルの推論値を利用した分析は「PyMC実装」章をご覧ください。

歯のきれいな猫のイラスト：「いらすとや」さんより

PyMC実装

Let's enjoy PyMC & Python !

準備・データ確認

１．インポート

### インポート

# ユーティリティ
import pickle

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

# PyMC
import pymc as pm
import bambi as bmb
import arviz as az

# 統計処理
import statsmodels.formula.api as smf
import lmdiag

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

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

２．データの読み込みと前処理
csvファイル「cvcps.csv」をpandasのデータフレームに読み込みます。

### データの読み込み
data = pd.read_csv('cvcps.csv')
display(data)

【実行結果】
項目内容の概要です。
・id：回答者ID
・qol：QOLの回答値（目的変数）
・v1～v5：CVCPS５項目の回答値

３．データの概観の確認
要約統計量を確認しましょう。

### 要約統計量の表示
data.describe().round(2)

【実行結果】
中央値はほぼ５です。回答の中間値は４なので若干高い値になっているようです。

６つの回答値を可視化しましょう。
ヴァイオリンプロットを描いてみます。

### 質問６項目のヴァイオリンプロットの描画
sns.violinplot(data.iloc[:, 1:], palette='Pastel1', linecolor='gray')
plt.grid(lw=0.5)

【実行結果】
目的変数 qol と形状がよく似ているのが v3 （治療選択）と v5 （健康情報意識）のようです。

変数間の相関をヒートマップで可視化しましょう。

### 質問６項目の相関係数・ヒートマップの描画
sns.heatmap(data.iloc[:, 1:].corr(), vmin=0, vmax=1, cmap='Greens', alpha=0.7,
            annot=True);

【実行結果】
目的変数 qol と最も高い相関関係があるのは v3 （相関係数 0.62）です。

モデルの構築

Formula 構文でPyMCモデルを構築できる「Bambi」を使うことにします。

モデルの数式表現
目指したいBambiのモデルの「なんちゃって数式」表記です。

$$
qol \sim v_1 + v_2 + v_3 + v_4 + v_5
$$

上記の Formula は暗黙的に切片 Intercept を含んでいます。

１．モデルの定義
数式表現を実直にモデル記述します。

### モデルの定義
model_bmb = bmb.Model(
        formula='qol ~ v1 + v2 + v3 + v4 + v5',     # フォーミュラ式
        data=data,                                  # データ
        family='gaussian',                          # （省略可）誤差分布
        link='identity',                            # （省略可）リンク関数
)

２．モデルの外観の確認

# モデルの表示
model_bmb

【実行結果】
Bambi は独自ロジックで事前分布を割当てしてくれます。
切片 Intercept および（偏）回帰係数$${v_1}$$～$${v_5}$$は正規分布に従う設定になっています。
また、目的変数 qol が従う正規分布の標準偏差 sigma は半ｔ分布に従う設定になりました。

### モデルの可視化
model_bmb.build()
model_bmb.graph()

【実行結果】
確率変数が横並びです！

３．事後分布からのサンプリング
Bambi モデルに対して .fit することでMCMCサンプリングを実行できます。
乱数生成数はテキストよりも少ないです。
nuts_sampler='numpyro'とすることで、numpyroをNUTSサンプラーに利用できます。
処理時間はおよそ40秒でした。

### 事後分布からのサンプリング ※NUTSサンプラーにnumpyroを使用 40秒
# テキスト: iter=2000, warmup=1000, chains=4, thin=10
idata_bmb = model_bmb.fit(draws=1000, tune=1000, chains=4, target_accept=0.8, 
                          nuts_sampler='numpyro', random_seed=1234)

【実行結果】省略

４．サンプリングデータの確認
$${\hat{R}}$$、トレースプロットを確認します。
事後分布の収束確認は$${\hat{R} \leq 1.1}$$とします。
ここでは$${\hat{R} > 1.01}$$のパラメータが無いことを確認します。

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

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

【実行結果】
$${\hat{R} > 1.01}$$のパラメータは０件でした。
全てのパラメータが$${\hat{R} \leq 1.1}$$であることを確認できました。

ざっくり事後分布サンプリングデータの要約統計量とトレースプロットを確認します。

### 推論データの要約統計情報の表示
pm.summary(idata_bmb, hdi_prob=0.95, round_to=3)

【実行結果】

トレースプロットで事後分布サンプリングデータの様子を確認します。

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

【実行結果】
左のグラフから、４本のマルコフ連鎖chainがほぼ同じ分布であることが分かります。
右のグラフでは線が満遍なく描画されています。
収束していると思われます。

５．推論データの保存
推論データを再利用する（かもしれない）場合に備えてファイルに保存しましょう。
idata_bmb をpickleで保存します。

### idataの保存 pickle
file = r'idata_bmb_ch14.pkl'
with open(file, 'wb') as f:
    pickle.dump(idata_bmb, f)

読み込みコードは次のとおりです。

### idataの読み込み pickle
file = r'idata_bmb_ch14.pkl'
with open(file, 'rb') as f:
    idata_bmb_load = pickle.load(f)

モデルの分析

パラメータの事後統計量を算出します。
テキストの図14.2に相当します。

### パラメータの推定結果の表示 ※表14.2に相当

# 事後統計量算出関数の定義
def calc_stats(x):
    # 95%信用区間の算出
    ci = np.quantile(x, q=[0.025, 0.975])
    # 戻り値：EAP、post.sd, 2.5%CI, 97.5%CI
    return np.mean(x), np.std(x), ci[0], ci[1]

## 初期値設定
# パラメータ名の設定
params = [f'v{i}' for i in range(1, 6)]
names = ['助言・指示の遵守', '積極的コミットメント', '治療選択', '私的治療等会話',
         '健康情報意識']
# 空のデータフレームの作成
stats_df1 = pd.DataFrame()

## 事後統計量の算出
# パラメータごとに事後統計量を算出してデータフレームに追加する処理を繰り返す
for p in params:
    # 事後統計量の算出
    tmp = pd.DataFrame({p: calc_stats(idata_bmb.posterior[p].data.flatten())})
    # データフレームの結合
    stats_df1 = pd.concat([stats_df1, tmp], axis=1)

## データフレームの最終化・表示：行列の転置.T → 列名の設定.set_axis
stats_df1 = stats_df1.T.set_axis(['事後平均値', '標準偏差', '2.5%CI', '97.5%CI'],
                                 axis=1)
display(stats_df1.round(2))

【実行結果】
テキストの統計値とほぼ同じ結果になりました。

パラメータの事後分布を可視化しましょう。
テキストの図14.1に相当します。

### パラメータの事後分布の描画 ※図14.1に相当

# 描画領域の指定（figレベル）
plt.figure(figsize=(8, 5))
# パラメータごとに描画処理を繰り返し処理
for i, param in enumerate(params):
    # パラメータの事後分布サンプリングデータの取り出し
    tmp = idata_bmb.posterior[param].data.flatten()
    # 描画領域の指定（axesレベル）
    ax = plt.subplot(2, 3, i+1)
    # ヒストグラム＋KDE曲線の描画
    hist = sns.histplot(tmp, kde=True, stat='density', bins=30, edgecolor='white',
                        ax=ax)
    # KDE最大値のxの取得
    kde = hist.get_lines()[0].get_data()
    kde_max_idx = np.argmax(kde[1])
    kde_max_x = kde[0][kde_max_idx]
    # 修飾：x=0の赤点線、タイトル
    ax.vlines(kde_max_x, 0, kde[1][kde_max_idx], color='tab:red', ls='--')
    ax.set(title=f'{param}: {names[i]}\n KDE最大値:{kde_max_x:.2f}')

# axes描画間隔の調整
plt.tight_layout()
plt.show()

【実行結果】

【分析～テキストにならって】
テキストの解説のとおり「QOLに最も大きな影響を与えている可能性があるのは v3（治療選択）であった。その次に v1（助言・指示の遵守）の効果が大きい」結果になりました。
患者自身がしっかりと考えて積極的に治療を「選ぶ」ことと、歯科医師からの助言・指示を守れることが、歯科通院による生活の質向上に影響するようです。
テキストのまとめのとおり、歯科医院は「患者が実行可能で理解できるような指導を行うことが必要」でしょうし、患者が治療に関心をもって積極的に治療に選択・参加できるよう、歯科医師による適度な情報開示が大切だな、って感じました。

以上で第14章は終わり、ません！
アディショナルタイムに突入です！

モデル２：PyMCでモデリング

せっかくですので PyMC のモデリングも楽しみましょう！

モデルの数式表現
目指したいPyMCのモデルの雰囲気を混ぜた「なんちゃって数式」表記です。

$$
\begin{align*}
\sigma &\sim \text{Uniform}\ (\text{lower}=0,\ \text{upper}=100)\\
\beta &\sim \text{Uniform}\ (\text{lower}=-100,\ \text{upper}=100,\ \text{dims}=xdims)\\
\mu &= X @ \beta \\
likelihood &\sim \text{Normal}\ (\text{mu}=\mu, \text{simga}=\sigma) \\
\end{align*}
$$

１．モデルの定義
数式表現を実直にモデル記述します。

### モデルの定義

### データの準備
# 説明変数：１列目に定数1を設定（切片用）、２列目以降はv1～v5
x_data = pd.concat([pd.Series(np.ones(len(data)), name='Intercept'),
                    data.iloc[:, 2:]], axis=1)
# 説明変数の名称リスト：定数＋v1～v5
params2 = ['Intercept'] + params

### モデルの定義
with pm.Model() as model2:
    ## coordの定義
    # データのインデックス
    model2.add_coord('data', values=data.index, mutable=True)
    # 説明変数
    model2.add_coord('xDims', values=params2, mutable=True)
    
    ## dataの定義
    # 目的変数：QOL
    y = pm.ConstantData('y', value=data['qol'].values, dims='data')
    # 説明変数：質問項目v1～v5
    X = pm.ConstantData('X', value=x_data.values, dims=('data', 'xDims'))
    
    ## 事前分布：一様分布
    sigma = pm.Uniform('sigma', lower=0, upper=100)
    beta = pm.Uniform('beta', lower=-100, upper=100, dims='xDims')
    
    ## muの算出 @は行列積（内積）
    mu = pm.Deterministic('mu', X @ beta, dims='data')
    
    ## 尤度：正規分布
    likelihood = pm.Normal('likelihood', mu=mu, sigma=sigma, observed=y,
                             dims='data')

【モデル注釈】

  • coordの定義
    座標に名前を付けたり、その座標が取りうる値を設定できます。
    今回は次の２つを設定しました。

    • データの行の座標：名前「data」、値「行インデックス」

    • 説明変数の座標：名前「xdims」、値「切片を含む説明変数の名称」

  • dataの定義
    次の４つを設定しました。

    • QOLの回答値：y （目的変数）

    • 切片を含む説明変数：X

  • パラメータの事前分布

    • モデルの数式表現のとおりです。

  • 尤度

    • 平均パラメータが重回帰式$${\mu=X@\beta}$$、標準偏差パラメータが$${\sigma}$$の正規分布です。

２．モデルの外観の確認

### モデルの表示
model2

【実行結果】
とてもシンプルです。

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

【実行結果】

３．事後分布からのサンプリング
乱数生成数はテキストよりも少ないです。
nuts_sampler='numpyro'とすることで、numpyroをNUTSサンプラーに利用できます。
処理時間はおよそ10秒でした。

### 事後分布からのサンプリング ※NUTSサンプラーにnumpyroを使用　10秒
# テキスト: iter=2000, warmup=1000, chains=4, thin=10
with model2:
    idata2 = pm.sample(draws=1000, tune=1000, chains=4, target_accept=0.8,
                       nuts_sampler='numpyro', random_seed=123)

【実行結果】省略

４．サンプリングデータの確認
$${\hat{R}}$$、トレースプロットを確認します。
事後分布の収束確認は$${\hat{R} \leq 1.1}$$とします。
ここでは$${\hat{R} > 1.01}$$のパラメータが無いことを確認します。

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

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

【実行結果】
$${\hat{R} > 1.01}$$のパラメータは０件でした。
全てのパラメータが$${\hat{R} \leq 1.1}$$であることを確認できました。

ざっくり事後分布サンプリングデータの要約統計量とトレースプロットを確認します。

### 推論データの要約統計情報の表示
pm.summary(idata2, hdi_prob=0.95, round_to=3)

【実行結果】

切片、（偏）回帰係数を見てみます。

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

【実行結果】

トレースプロットで事後分布サンプリングデータの様子を確認します。

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

【実行結果】
左のグラフから、４本のマルコフ連鎖chainがほぼ同じ分布であることが分かります。
右のグラフでは線が満遍なく描画されています。
収束していると思われます。

５．推論データの保存
推論データを再利用する（かもしれない）場合に備えてファイルに保存しましょう。
idata2 をpickleで保存します。

### idataの保存 pickle
file = r'idata2_ch14.pkl'
with open(file, 'wb') as f:
    pickle.dump(idata2, f)

読み込みコードは次のとおりです。

### idataの読み込み pickle
file = r'idata2_ch14.pkl'
with open(file, 'rb') as f:
    idata2_load = pickle.load(f)

モデル２の分析

アウトプットを掲載するだけに留めます。
パラメータの事後統計量を算出します。
テキストの図14.2に相当します。

### パラメータの推定結果の表示 ※表14.2に相当

# 空のデータフレームの作成
stats_df2 = pd.DataFrame()

## 事後統計量の算出
# パラメータごとに事後統計量を算出してデータフレームに追加する処理を繰り返す
for i, p in enumerate(params):
    # 事後統計量の算出
    idata_tmp = idata2.posterior.beta.stack(sample=('chain', 'draw')).data[i+1]
    tmp_df = pd.DataFrame({p: calc_stats(idata_tmp)})
    # データフレームの結合
    stats_df2 = pd.concat([stats_df2, tmp_df], axis=1)

## データフレームの最終化・表示：行列の転置.T → 列名の設定.set_axis
stats_df2 = stats_df2.T.set_axis(['事後平均値', '標準偏差', '2.5%CI', '97.5%CI'],
                                 axis=1)
display(stats_df2.round(2))

【実行結果】
Bambiモデルとほぼ同じ結果です。

（参考：Bambiモデルの事後統計量）

パラメータの事後分布を可視化しましょう。
テキストの図14.1に相当します。

### パラメータの事後分布の描画 ※図14.1に相当

# 描画領域の指定（figレベル）
plt.figure(figsize=(8, 5))
# パラメータごとに描画処理を繰り返し処理
for i, param in enumerate(params):
    # パラメータの事後分布サンプリングデータの取り出し
    tmp = idata2.posterior.beta.stack(sample=('chain', 'draw')).data[i+1]
    # 描画領域の指定（axesレベル）
    ax = plt.subplot(2, 3, i+1)
    # ヒストグラム＋KDE曲線の描画
    hist = sns.histplot(tmp, kde=True, stat='density', bins=30, edgecolor='white',
                        ax=ax)
    # KDE最大値のxの取得
    kde = hist.get_lines()[0].get_data()
    kde_max_idx = np.argmax(kde[1])
    kde_max_x = kde[0][kde_max_idx]
    # 修飾：x=0の赤点線、タイトル
    ax.vlines(kde_max_x, 0, kde[1][kde_max_idx], color='tab:red', ls='--')
    ax.set(title=f'{param}: {names[i]}\n KDE最大値:{kde_max_x:.2f}')

# axes描画間隔の調整
plt.tight_layout()
plt.show()

【実行結果】

おまけ

statsmodels で統計的な重回帰分析を行ってみます。

### statsmodelsのformula APIで線形回帰を実施
# モデルの定義・学習の実行
lm = smf.ols(formula='qol ~ v1 + v2 + v3 + v4 + v5', data=data).fit()
# 実行結果の表示
print(lm.summary())

【実行結果】
（偏）回帰係数 coef は概ねベイズ推論の平均値に近いです。

（偏）回帰係数と95%信頼区間（ベイズ信用区間では無いです）を表示します。

### 回帰係数と95%信頼区間の表示
lm_summary = lm.summary2().tables[1]
display(lm_summary.round(2))

【実行結果】
Bambiモデルの結果と値がよく似ています。

（参考：Bambiモデルの事後統計量）

可視化にもチャレンジします。

### 回帰係数・標準誤差を用いた正規分布の確率密度関数の描画

# 初期値設定：x軸目盛ラベル
x = np.linspace(-0.2, 0.8, 1001)
# 描画領域の指定（figレベル）
plt.figure(figsize=(8, 5))
for i in range(1, 6):
    # 平均・標準偏差（標準誤差）の取得
    loc, scale = lm_summary.iloc[i, :2]
    # 正規分布の確率密度関数の取得
    density = stats.norm.pdf(x=x, loc=loc, scale=scale)
    # 描画領域の指定（axesレベル）
    ax = plt.subplot(2, 3, i)
    # 正規分布の確率密度関数の描画
    ax.plot(x, density)
    # 最大値のxの取得
    max_idx = np.argmax(density)
    max_x = x[max_idx]
    # 修飾：x=0の赤点線、タイトル
    ax.vlines(max_x, 0, density[max_idx], color='tab:red', ls='--')
    ax.set(title=f'{params2[i]}: {names[i-1]}\n 確率密度最大値:{max_x:.2f}',
           xlim=(loc-scale*4, loc+scale*4))
# axes描画間隔の調整
plt.tight_layout()
plt.show()

【実行結果】
鏡写しのようなフォルムになりました・・・

残差の様子を確認しましょう
lmdiag ライブラリを利用します。

### 残差プロットの描画
plt.figure(figsize=(7, 7))
lmdiag.plot(lm);

【実行結果】

以上で第14章は終了です。

おわりに

QOL

以前からQOLが「ほぼ」医療・福祉の分野の言葉として扱われていることに違和感を覚えていました。
生活の質は人生の様々なシーンで大切にしたいことであって、医療・福祉に特化する概念では無いんじゃないかって。

ところで、ベイズでQOLは向上するでしょうか？
広く一般的に向上するかどうかは分からないです。
ただ個人的には、ベイズ推論や機械学習などを学ぶ過程で、自己肯定感が高くなっているように感じています。
自己肯定感が高まるとQOLも向上するそうなので、私はベイズでQOLが向上しているのかもしれません。

おわり

シリーズの記事

次の記事

前の記事

目次


ブログの紹介

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

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

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

新生活をたのしく

44,955件

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