StanとRでベイズ統計モデリングをPyMC Ver.5で写経~第5章「5.1 重回帰」
第5章「基本的な回帰とモデルのチェック」
書籍の著者 松浦健太郎 先生
この記事は、テキスト第5章「基本的な回帰とモデルのチェック」の5.1節「重回帰」の PyMC5写経 を取り扱います。
今回は特にグラフを描画するコードがたくさん出現します。
Python化するのにとても時間がかかりました(汗)
はじめに
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を動かすまでの準備」章をご覧ください。
5.1 重回帰
インポート
### インポート
# 数値・確率計算
import pandas as pd
import numpy as np
import scipy.stats as stats
# PyMC
import pymc as pm
import arviz as az
# 描画
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import seaborn as sns
plt.rcParams['font.family'] = 'Meiryo'
# ワーニング表示の抑制
import warnings
warnings.simplefilter('ignore')
サンプルコードのデータを読み込みます。
### データの読み込み ◆データファイル5.1 data-attendance-1.txtの構成
# A:バイト好き区分(1:好き), Score:学問の興味の強さ(0~200), Y:授業の出席率(1年間)
data = pd.read_csv('./data/data-attendance-1.txt')
print('data.shape: ', data.shape)
display(data.head())
【実行結果】
5.1.2 データの分布の確認
散布図行列を描画します。図5.1相当です。
なるべくテキストに寄せるように努力しました。
### 散布図行列の描画 ◆図5.1
# 凡例非表示・・・描画関数の引数に legend=None を追加する
## 描画領域の指定
fig, ax = plt.subplots(3, 3, figsize=(10, 10))
ax = ax.ravel() # 1次元でaxesを指定したいので
## 番地0,0:ヒストグラムの描画(棒グラフを使用)
bar_A = data.A.value_counts().sort_index()
sns.barplot(ax=ax[0], x=bar_A.index, y=bar_A, hue=bar_A.index, palette='tab10',
alpha=0.5, ec='white')
ax[0].set(ylabel='A', xlabel=None)
ax[0].grid(lw=0.5)
## 番地0,1:スピアマンの順位相関係数の描画
ax[1].set_axis_off()
corr1, pval1 = stats.spearmanr(data.Score, data.A)
ax[1].text(x=0.4, y=0.4, s=round(corr1 * 100), fontsize=50)
## 番地0,2:スピアマンの順位相関係数の描画
ax[2].set_axis_off()
corr2, pval2 = stats.spearmanr(data.Y, data.A)
ax[2].text(x=0.2, y=0.4, s=round(corr2 * 100), fontsize=50)
ax[2].text(x=0, y=0.9, s='スピアマンの順位相関係数', fontsize=14)
## 番地1,0:箱ひげ図+スウォームプロットの描画
sns.boxplot(ax=ax[3], x=data.A, y=data.Score, hue=data.A, fill=False,
legend=None)
sns.swarmplot(ax=ax[3], x=data.A, y=data.Score, hue=data.A, size=8, alpha=0.5)
ax[3].set(xlabel=None)
ax[3].grid(lw=0.5)
## 番地1,1:ヒストグラムの描画
sns.histplot(ax=ax[4], data=data, x='Score', hue='A', bins=10, kde=True,
ec='white')
ax[4].set(xlabel=None, ylabel=None)
ax[4].grid(lw=0.5)
## 番地1,2:スピアマンの順位相関係数の描画
ax[5].set_axis_off()
corr3, pval3 = stats.spearmanr(data.Y, data.Score)
ax[5].text(x=0.3, y=0.4, s=round(corr3 * 100), fontsize=50)
## 番地2,0:箱ひげ図+スウォームプロットの描画
sns.boxplot(ax=ax[6], x=data.A, y=data.Y, hue=data.A, fill=False,
legend=None)
sns.swarmplot(ax=ax[6], x=data.A, y=data.Y, hue=data.A, size=8, alpha=0.5)
ax[6].grid(lw=0.5)
## 番地2,1:散布図の描画
sns.scatterplot(ax=ax[7], data=data, x='Score', y='Y', hue='A', size='A',
sizes=(80, 80), alpha=0.5)
ax[7].set(ylabel=None)
ax[7].grid(lw=0.5)
## 番地2,2:ヒストグラムの描画
sns.histplot(ax=ax[8], data=data, x='Y', hue='A', bins=10, kde=True, ec='white')
ax[8].set(ylabel=None)
ax[8].grid(lw=0.5)
plt.tight_layout();
【実行結果】
5.1.5 Stanで実装
PyMC Ver.5 で実装します。
モデルの定義です。
### モデルの定義 ◆model5-3.stan
with pm.Model() as model:
### データ関連定義
## coordの定義
model.add_coord('data', values=data.index, mutable=True)
## dataの定義
# 目的変数 Y
Y = pm.ConstantData('Y', value=data['Y'].values, dims='data')
# 説明変数 A
A = pm.ConstantData('A', value=data['A'].values, dims='data')
# 説明変数 Score / 200
Score = pm.ConstantData('Score', value=data['Score'].values / 200,
dims='data')
### 事前分布
b1 = pm.Uniform('b1', lower=-10, upper=10)
b2 = pm.Uniform('b2', lower=-10, upper=10)
b3 = pm.Uniform('b3', lower=-10, upper=10)
sigma = pm.Uniform('sigma', lower=0, upper=10)
### 線形予測子
mu = pm.Deterministic('mu', b1 + b2 * A + b3 * Score, dims='data')
### 尤度関数
obs = pm.Normal('obs', mu=mu, sigma=sigma, observed=Y, dims='data')
### 計算値
yPred = pm.Normal('yPred', mu=mu, sigma=sigma, dims='data')
モデルの定義内容を見ます。
### モデルの表示
model
【実行結果】
### モデルの可視化
pm.model_to_graphviz(model)
【実行結果】
5.1.6 データのスケーリング
PythonでMCMCを実行します。
### 事後分布からのサンプリング 20秒 ◆run-model5-3.R
with model:
idata = pm.sample(draws=1000, tune=1000, chains=4, target_accept=0.8,
nuts_sampler='numpyro', random_seed=1234)
【実行結果】省略
5.1.7 推定結果の解釈
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 = ['b1', 'b2', 'b3', 'sigma', 'mu']
pm.summary(idata, hdi_prob=0.95, var_names=var_names, round_to=3)
【実行結果】
テキスト59ページ掲載のベイズ信頼区間(50%, 95%)を含む事後統計量を表示します。
### パラメータの要約を確認 ◆テキスト59ページの要約統計量
# mean,sd,2.5%,25%,50%,75%,97.5%パーセンタイル点をデータフレーム化する関数の定義
def make_stats_df(y):
probs = [2.5, 25, 50, 75, 97.5]
columns = ['mean', 'sd'] + [str(s) + '%' for s in probs]
quantiles = pd.DataFrame(np.percentile(y, probs, axis=0).T, index=y.columns)
tmp_df = pd.concat([y.mean(axis=0), y.std(axis=0), quantiles], axis=1)
tmp_df.columns=columns
return tmp_df
# 要約統計量の算出・表示
vars = ['b1', 'b2', 'b3', 'sigma']
param_samples = idata.posterior[vars].to_dataframe().reset_index(drop=True)
display(make_stats_df(param_samples).round(2))
【実行結果】
トレースプロットを描画します。
### トレースプロットの表示
var_names = ['b1', 'b2', 'b3', 'sigma', 'mu', 'yPred']
pm.plot_trace(idata, compact=True, var_names=var_names)
plt.tight_layout();
【実行結果】
5.1.8 図によるモデルのチェック
■ 図5.2
図5.2のアルバイト好き区分別のベイズ予測区間を描画します。
出席率$${Y}$$の予測値は、MCMCサンプリングデータを用いず、このコード内で x 軸のScoreの範囲$${[50, 200]}$$に対応させて計算しています。
### 説明変数Scoreと目的変数の予測分布yPredをアルバイト好き区分ごとに描画 ◆図5.2
## 描画用データの作成
# 乱数生成器の設定
rng = np.random.default_rng(seed=123)
# X軸のScoreについて50~200の連続値を1001個作成
x_scores = np.linspace(50, 200, 1001)
# MCMCサンプリングデータからb1~sigmaを取り出し
b1_samples = idata.posterior.b1.stack(sample=('chain', 'draw')).data
b2_samples = idata.posterior.b2.stack(sample=('chain', 'draw')).data
b3_samples = idata.posterior.b3.stack(sample=('chain', 'draw')).data
sigma_samples = idata.posterior.sigma.stack(sample=('chain', 'draw')).data
# アルバイト0とアルバイト1ごとに平均値muを計算(回帰式)
mu_A0 = np.array([b1_samples + b2_samples * 0 + b3_samples * x_score / 200
for x_score in x_scores])
mu_A1 = np.array([b1_samples + b2_samples * 1 + b3_samples * x_score / 200
for x_score in x_scores])
# アルバイト0とアルバイト1ごとに正規分布乱数でYの予測分布を計算
y_pred_A0 = rng.normal(loc=mu_A0, scale=sigma_samples)
y_pred_A1 = rng.normal(loc=mu_A1, scale=sigma_samples)
# アルバイト0とアルバイト1ごとに予測分布のYの80%CIを計算
y_pred_A0_80ci = np.quantile(y_pred_A0, q=[0.1, 0.9], axis=1)
y_pred_A1_80ci = np.quantile(y_pred_A1, q=[0.1, 0.9], axis=1)
## 描画処理
# 観測値の散布図の描画
sns.scatterplot(data=data, x='Score', y='Y', hue='A', size='A', sizes=(80, 80),
alpha=0.5)
# アルバイト0のyの予測分布の中央値を描画
plt.plot(x_scores, np.median(y_pred_A0, axis=1), color='tab:blue')
# アルバイト0のyの予測分布の80%CIを描画
plt.fill_between(x_scores, y_pred_A0_80ci[0], y_pred_A0_80ci[1],
color='tab:blue', alpha=0.2)
# アルバイト1のyの予測分布の中央値を描画
plt.plot(x_scores, np.median(y_pred_A1, axis=1), color='tab:orange')
# アルバイト1のyの予測分布の80%CIを描画
plt.fill_between(x_scores, y_pred_A1_80ci[0], y_pred_A1_80ci[1],
color='tab:orange', alpha=0.2)
# 修飾
plt.title('出席率$Y$のベイズ80%予測区間')
plt.legend(bbox_to_anchor=(1, 1), title='アルバイト別\n$Y$の觀測値')
plt.grid(lw=0.5);
【実行結果】
■ 図5.3
図5.3の観測値と予測値のプロットを描画します。
matplotlib.pyplot の errorbar プロットを用いて、y軸方向に伸びるerrorbarを描画しました。
### 観測値と予測値のプロット ◆図5.3
## 描画用データの作成 yPredの個人別の中央値と80%区間を算出
# MCMCサンプリングデータからyPredを取り出し
y_pred_samples = idata.posterior.yPred.stack(sample=('chain', 'draw')).data
# サンプリングデータの10%,50%,90%パーセンタイル点を算出してデータフレーム化
y_pred_df = pd.DataFrame(
np.quantile(y_pred_samples, q=[0.1, 0.5, 0.9], axis=1).T,
columns=['10%', 'median', '90%'])
y_pred_df = pd.concat([data, y_pred_df], axis=1)
# 中央値と10%点の差、90%点と中央値の差を算出: errorbarで利用
y_pred_df['err_lower'] = y_pred_df['median'] - y_pred_df['10%']
y_pred_df['err_upper'] = y_pred_df['90%'] - y_pred_df['median']
# アルバイト0とアルバイト1に分離
y_pred_A0 = y_pred_df[y_pred_df['A']==0]
y_pred_A1 = y_pred_df[y_pred_df['A']==1]
## 描画処理
# 描画領域の指定
plt.figure(figsize=(6, 6))
ax = plt.subplot()
# アルバイト0の描画(エラーバー付き散布図)
ax.errorbar(y_pred_A0['Y'], y_pred_A0['median'],
yerr=[y_pred_A0['err_lower'], y_pred_A0['err_upper']],
color='tab:blue', alpha=0.7, marker='o', ms=10, linestyle='none',
label='0')
# アルバイト1の描画(エラーバー付き散布図)
ax.errorbar(y_pred_A1['Y'], y_pred_A1['median'],
yerr=[y_pred_A1['err_lower'], y_pred_A1['err_upper']],
color='tab:orange', alpha=0.7, marker='^', ms=10, linestyle='none',
label='1')
# 赤い対角線の描画
ax.plot([0, 0.5], [0, 0.5], color='red', ls='--')
# 修飾
ax.set(xlabel='Observed: $Y$の観測値', ylabel='Predicted: $Y$の予測値(中央値)',
title='$Y$ の観測値と予測値(中央値)のプロット 80%区間バー付き')
ax.legend(title='アルバイト好き', bbox_to_anchor=(1, 1))
ax.grid(lw=0.5);
【実行結果】
■ 図5.4左
図5.4左の個人ごとの誤差$${\varepsilon[n]}$$の分布を描画します。
seaborn の kdeplot を用いて分布を描画するとともに、kde_plot の線から KDEの値を取得して(get_lines()[i].get_data())、MAP推定値の取得に利用しています。
### 個人ごとのεの分布の描画 ◆図5.4左
## 描画用データの作成 yPredの個人別の中央値と80%区間を算出
# MCMCサンプリングデータからyPredを取り出し
mu_samples = idata.posterior.mu.stack(sample=('chain', 'draw')).data
# εの算出: 観測値Y - muの事後分布サンプリングデータ
epsilon_samples = np.array([data.loc[i, 'Y'] - mu_sample
for i, mu_sample in enumerate(mu_samples)])
## 図5.4右の描画に利用するMAP推定値リストの初期化
map_plot = []
## 描画処理
# 描画領域の指定
plt.figure(figsize=(15, 5))
# 個人ごとのε(ε[n])のKDEプロットの描画
for i, epsilon in enumerate(epsilon_samples):
# 色の設定
color = cm.tab20(i % 20)
# x軸の値を設定
x_ticks = np.linspace(epsilon.min(), epsilon.max(), 1001)
# KDEプロットの描画
kde_plot = sns.kdeplot(epsilon, color=color)
# KDEプロットで描画した線から、x軸,y軸の値を取得
kde_data = kde_plot.get_lines()[i].get_data()
# KDEプロットのyの最大値のインデックスを取得
kde_max_idx = kde_data[1].argmax()
# KDEプロットのyの最大値≒MAP推定値の垂直線の描画
plt.vlines(x=kde_data[0][kde_max_idx], ymin=0, ymax=kde_data[1][kde_max_idx],
color=color, ls='--')
# KDEプロットのyの最大値≒MAP推定値のy=0付近のx点の描画
plt.scatter((kde_data[0][kde_max_idx]), (-1), marker='s', s=20, color=color)
# 図5.4右用のMAP推定値の格納
map_plot.append(kde_data[0][kde_max_idx])
# 修飾
plt.xlabel('value')
plt.title('個人ごとの誤差 $\epsilon[n]$ の分布:KDEプロット')
plt.ylim(-3, 45)
plt.grid(lw=0.5)
【実行結果】
テキストには無いですが、$${\sigma}$$のKDEプロットを描画します。
$${\sigma}$$のMAP推定値算出の目的で、KDEプロットの線の値を取得したいからです。
### σのKDEプロット
sigma_plot = sns.kdeplot(sigma_samples)
sigma_kde_data = sigma_plot.get_lines()[0].get_data()
sigma_map = sigma_kde_data[0][sigma_kde_data[1].argmax()]
plt.vlines(sigma_map, 0, sigma_kde_data[1].max(), colors='tab:red', ls='--')
plt.title(f'標準偏差のMAP推定値: {sigma_map:.3f}');
【実行結果】
■ 図5.4右
図5.4右の$${\varepsilon[n]}$$のMAP推定値の分布を描画します。
### εのMAP推定値のヒストグラム・KDEプロット、σのMAP推定値による正規分布 ◆図5.4右
## 描画用データの作成
# 正規分布のx軸の値を設定
xvals = np.linspace(-0.16, 0.16, 1001)
# 正規分布の確率密度関数の算出
yvals = stats.norm.pdf(xvals, loc=0, scale=sigma_map)
## 描画処理
# σのMAP推定値を用いた正規分布の確率密度関数の描画
plt.plot(xvals, yvals, color='tab:orange', label='正規分布(0, $\sigma$)')
# εのMAP推定値のKDEプロットの描画
sns.kdeplot(np.array(map_plot), color='tab:blue',
label='$\epsilon$のMAP推定値のKDE')
# εのMAP推定値のヒストグラムの描画
sns.histplot(map_plot, bins=15, color='tab:blue', ec='white', alpha=0.7,
label='$\epsilon$のMAP推定値')
# 修飾
plt.title('$\epsilon[n]$ のMAP推定値のヒストグラム・KDE,\n'
'$\sigma$ のMAP推定値を用いた正規分布の確率密度関数')
plt.xlabel('value')
plt.legend(bbox_to_anchor=(1.45, 1))
plt.grid(lw=0.5);
【実行結果】
■ 図5.5
図5.5のMCMCサンプルの散布図行列を描画します。
seaborn の pairplot を利用しており、テキストの描画内容よりも簡略化しています。
### MCMCサンプルの散布図行列の描画 ◆図5.5を簡素化
## 描画用データの作成
# MCMCサンプリングデータからmu1, mu50を取り出し
mu1_samples = idata.posterior.mu.stack(sample=('chain', 'draw'))[0]
mu50_samples = idata.posterior.mu.stack(sample=('chain', 'draw'))[-1]
# 描画対象パラメータをデータフレーム化
plot_df = pd.DataFrame({
'b1': b1_samples, 'b2': b2_samples, 'b3' : b3_samples,
'sigma': sigma_samples, 'mu1': mu1_samples, 'mu50': mu50_samples})
## 描画処理
# 相関行列プロットの描画
g = sns.pairplot(plot_df, diag_kws={'kde': True, 'ec': 'white'})
# スピアマンの順位相関係数の表示のためのaxフラット化
ax = g.axes.ravel()
## スピアマンの順位相関係数の表示
# b1, b2
corr1, pval1 = stats.spearmanr(plot_df.b1, plot_df.b2)
ax[1].text(x=-0.11, y=0.24, s=round(corr1 * 100), fontsize=20, color='red')
# b1, b3
corr1, pval1 = stats.spearmanr(plot_df.b1, plot_df.b3)
ax[2].text(x=0.4, y=0.24, s=round(corr1 * 100), fontsize=20, color='red')
# b1, sigma
corr1, pval1 = stats.spearmanr(plot_df.b1, plot_df.sigma)
ax[3].text(x=0.078, y=0.24, s=round(corr1 * 100), fontsize=20, color='red')
# b1, mu1
corr1, pval1 = stats.spearmanr(plot_df.b1, plot_df.mu1)
ax[4].text(x=0.3, y=0.24, s=round(corr1 * 100), fontsize=20, color='red')
# b1, mu50
corr1, pval1 = stats.spearmanr(plot_df.b1, plot_df.mu50)
ax[5].text(x=0.18, y=0.24, s=round(corr1 * 100), fontsize=20, color='red')
# b2, b3
corr1, pval1 = stats.spearmanr(plot_df.b2, plot_df.b3)
ax[8].text(x=0.45, y=-0.1, s=round(corr1 * 100), fontsize=20, color='red')
# b2, sigma
corr1, pval1 = stats.spearmanr(plot_df.b2, plot_df.sigma)
ax[9].text(x=0.078, y=-0.1, s=round(corr1 * 100), fontsize=20, color='red')
# b2, mu1
corr1, pval1 = stats.spearmanr(plot_df.b2, plot_df.mu1)
ax[10].text(x=0.29, y=-0.1, s=round(corr1 * 100), fontsize=20, color='red')
# b2, mu50
corr1, pval1 = stats.spearmanr(plot_df.b2, plot_df.mu50)
ax[11].text(x=0.18, y=-0.1, s=round(corr1 * 100), fontsize=20, color='red')
# b3, sigma
corr1, pval1 = stats.spearmanr(plot_df.b3, plot_df.sigma)
ax[15].text(x=0.078, y=0.45, s=round(corr1 * 100), fontsize=20, color='red')
# b3, mu1
corr1, pval1 = stats.spearmanr(plot_df.b3, plot_df.mu1)
ax[16].text(x=0.29, y=0.45, s=round(corr1 * 100), fontsize=20, color='red')
# b3, mu50
corr1, pval1 = stats.spearmanr(plot_df.b3, plot_df.mu50)
ax[17].text(x=0.175, y=0.45, s=round(corr1 * 100), fontsize=20, color='red')
# sigma, mu1
corr1, pval1 = stats.spearmanr(plot_df.sigma, plot_df.mu1)
ax[22].text(x=0.31, y=0.075, s=round(corr1 * 100), fontsize=20, color='red')
# sigma, mu50
corr1, pval1 = stats.spearmanr(plot_df.sigma, plot_df.mu50)
ax[23].text(x=0.19, y=0.075, s=round(corr1 * 100), fontsize=20, color='red')
# mu1, mu50
corr1, pval1 = stats.spearmanr(plot_df.mu1, plot_df.mu50)
ax[29].text(x=0.18, y=0.29, s=round(corr1 * 100), fontsize=20, color='red');
【実行結果】
5.1 節は以上です。
シリーズの記事
次の記事
前の記事
目次
ブログの紹介
note で7つのシリーズ記事を書いています。
ぜひ覗いていってくださいね!
1.のんびり統計
統計検定2級の問題集を手がかりにして、確率・統計をざっくり掘り下げるブログです。
雑談感覚で大丈夫です。ぜひ覗いていってくださいね。
統計検定2級公式問題集CBT対応版に対応しています。
Python、EXCELのサンプルコードの配布もあります。
2.実験!たのしいベイズモデリング1&2をPyMC Ver.5で
書籍「たのしいベイズモデリング」・「たのしいベイズモデリング2」の心理学研究に用いられたベイズモデルを PyMC Ver.5で描いて分析します。
この書籍をはじめ、多くのベイズモデルは R言語+Stanで書かれています。
PyMCの可能性を探り出し、手軽にベイズモデリングを実践できるように努めます。
身近なテーマ、イメージしやすいテーマですので、ぜひぜひPyMCで動かして、一緒に楽しみましょう!
3.実験!岩波データサイエンス1のベイズモデリングをPyMC Ver.5で
書籍「実験!岩波データサイエンスvol.1」の4人のベイジアンによるベイズモデルを PyMC Ver.5で描いて分析します。
この書籍はベイズプログラミングのイロハをざっくりと学ぶことができる良書です。
楽しくPyMCモデルを動かして、ベイズと仲良しになれた気がします。
みなさんもぜひぜひPyMCで動かして、一緒に遊んで学びましょう!
4.楽しい写経 ベイズ・Python等
ベイズ、Python、その他の「書籍の写経活動」の成果をブログにします。
主にPythonへの翻訳に取り組んでいます。
写経に取り組むお仲間さんのサンプルコードになれば幸いです🍀
5.RとStanではじめる心理学のための時系列分析入門 を PythonとPyMC Ver.5 で
書籍「RとStanではじめる心理学のための時系列分析入門」の時系列分析をPythonとPyMC Ver.5 で実践します。
この書籍には時系列分析のテーマが盛りだくさん!
時系列分析の懐の深さを実感いたしました。
大好きなPythonで楽しく時系列分析を学びます。
6.データサイエンスっぽいことを綴る
統計、データ分析、AI、機械学習、Pythonのコラムを不定期に綴っています。
統計・データサイエンス書籍にまつわる記事が多いです。
「統計」「Python」「数学とPython」「R」のシリーズが生まれています。
7.Python機械学習プログラミング実践記
書籍「Python機械学習プログラミング PyTorch & scikit-learn編」を学んだときのさまざまな思いを記事にしました。
この書籍は、scikit-learnとPyTorchの教科書です。
よかったらぜひ、お試しくださいませ。
最後までお読みいただきまして、ありがとうございました。