PythonでStanを使用する方法

プログラミング
2024.08.09
記事内に広告が含まれています。
スポンサーリンク

はじめに

brian
brian

ベイズモデリングって何だろう？PythonとStanを使ってその基本から学べるガイドを用意しました！統計やモデリングが初めての方でも、少しずつ理解できるよう丁寧に解説しているので、ぜひ気軽に読んでみてくださいね！

Stanとは何か？

Stanは、ベイズ統計モデリングと機械学習のためのプログラミング言語およびソフトウェアライブラリです。多くの現実世界の問題は複雑で不確実性が伴うため、データに対する柔軟で強力なモデリング手法が必要です。Stanは、確率モデルを定義し、ベイズ推定やマルコフ連鎖モンテカルロ（MCMC）サンプリングを通じてこれらのモデルを実行するためのツールを提供します。

Stanは、R、Python、Julia、Matlabなど、複数のプログラミング言語からアクセス可能です。その中でもPythonは、科学計算とデータ分析の強力なエコシステムを持っており、Stanを利用する際に非常に便利です。

PythonでStanを使う利点

PythonでStanを使用することにはいくつかの利点があります。まず、Pythonは幅広いライブラリとツールが揃っており、データの前処理や結果の可視化が容易です。また、Pythonのコミュニティは活発で、多くのドキュメントやサポートが利用可能です。

StanのPythonラッパーであるpystanやcmdstanpyを利用することで、Stanの機能をフルに活用しながら、Pythonの利便性を享受できます。これにより、複雑な統計モデリングがよりシームレスに行えます。

ベイズ推定とStanの役割

ベイズ推定は、データに基づいて不確実性を考慮しながらモデルを構築する手法です。従来の統計手法では、モデルのパラメータを一点として推定しますが、ベイズ推定ではパラメータに対する確率分布を求めます。これにより、データに対するより柔軟な解釈が可能になります。

Stanは、ベイズ推定を実行するための強力なツールであり、複雑な確率モデルを定義し、効率的に計算することができます。Stanを使えば、複雑な階層モデルや非線形モデルなど、従来の手法では扱いにくい問題にも対応できるのです。

スポンサーリンク

PythonでStanをセットアップする

必要なライブラリのインストール

PythonでStanを使用するためには、いくつかの専用ライブラリをインストールする必要があります。代表的なものとして、pystanとcmdstanpyが挙げられます。これらのライブラリは、PythonからStanを呼び出し、モデルのコンパイルやサンプリングを行うためのインターフェースを提供します。

  • PyStan: PyStanはStanのPythonラッパーであり、StanのモデルをPythonから直接操作することができます。インストールは次のコマンドで行えます。
Bash
pip install pystan
  • CmdStanPy: CmdStanPyは、Stanのコマンドラインインターフェース（CmdStan）をPythonから操作するためのライブラリです。CmdStanPyを使うことで、最新のStanの機能を活用できます。インストールは次のコマンドで行えます。
Bash
pip install cmdstanpy

インストール手順と注意点

ライブラリのインストールに際して、以下の点に注意する必要があります。

  1. Pythonのバージョン
    PyStanやCmdStanPyを使用するためには、Python 3.7以降が推奨されます。特に最新のStanの機能を使いたい場合は、Pythonのバージョンを確認しておくことが重要です。
  2. コンパイラの設定
    StanはC++で書かれているため、モデルのコンパイルにはC++コンパイラが必要です。通常、CmdStanPyを使用する場合には自動的にC++コンパイラがインストールされますが、システムに合わせて適切なコンパイラが設定されていることを確認する必要があります。
    • Windows: Windowsでは、Microsoft Visual StudioのBuild Toolsが必要です。
    • macOS: Xcode Command Line Toolsをインストールします。
    • Linux: g++やclangなどの標準的なC++コンパイラが必要です。
  3. CmdStanのセットアップ
    CmdStanPyを利用する場合、CmdStan自体もインストールする必要があります。CmdStanはスタンドアロンのプログラムで、Pythonからの操作にはこのプログラムが必須です。CmdStanPyはインストール時にCmdStanも自動的にダウンロードしますが、手動でのセットアップも可能です。
Bash
cmdstanpy.install_cmdstan()

サンプルプロジェクトの設定

ライブラリのインストールが完了したら、次にサンプルプロジェクトを設定してみましょう。ここでは、最もシンプルなモデルを使用して、Stanが正しくインストールされ、動作していることを確認します。

  1. モデルファイルの作成
    まず、simple_model.stanというファイルを作成し、以下のStanコードを記述します。
stan
data {
  int<lower=0> N;
  real y[N];
}
parameters {
  real mu;
}
model {
  mu ~ normal(0, 1);
  y ~ normal(mu, 1);
}

このモデルは、平均値muを推定するためのシンプルな正規分布モデルです。

  1. Pythonスクリプトの作成
    次に、このStanモデルをPythonで実行するためのスクリプトを作成します。以下は、pystanを使った例です。
Python
import pystan

model_code = """
data {
  int<lower=0> N;
  real y[N];
}
parameters {
  real mu;
}
model {
  mu ~ normal(0, 1);
  y ~ normal(mu, 1);
}
"""
stan_model = pystan.StanModel(model_code=model_code)

data = {'N': 10, 'y': [2.9, 3.1, 2.5, 2.8, 3.0, 3.2, 2.7, 3.0, 3.1, 2.9]}
fit = stan_model.sampling(data=data, iter=1000, chains=4)

print(fit)

このスクリプトを実行することで、サンプルデータに基づいてmuの値が推定され、結果が表示されます。

  1. 結果の確認
    スクリプトを実行し、結果が正しく出力されることを確認します。出力には、muの推定値や、その事後分布の要約が含まれます。
スポンサーリンク

Stanモデルの基本構造

Stanでは、確率モデルを定義するために独自のプログラミング言語を使用します。Stanのプログラムは、いくつかの主要なブロックで構成されており、これらのブロックを組み合わせることで複雑なモデルを定義します。この章では、Stanプログラムの基本構造を説明し、PythonからStanモデルを定義する方法を紹介します。

Stanプログラムの基本的な構成要素

Stanプログラムは、通常以下の4つの主要なブロックで構成されています。

  1. データブロック (data block)
    モデルに入力されるデータを定義するブロックです。ここで宣言される変数は、外部から与えられるデータを受け取ります。変数の型やサイズを指定し、モデル内で使用される前にデータの範囲や制約を明確にします。
stan
data {
  int<lower=0> N; // データポイントの数
  real y[N]; // 観測されたデータ
}
  1. パラメータブロック (parameters block)
    モデルが推定しようとする未知のパラメータを定義するブロックです。ここで宣言された変数は、ベイズ推定によってサンプルされ、結果として得られます。パラメータには範囲や制約を設定することも可能です。
stan
parameters {
  real mu; // 推定する平均値
}
  1. モデルブロック (model block)
    モデルの構造とデータに対する確率分布を定義するブロックです。このブロックでは、データがどのように生成されたかを記述します。具体的には、パラメータに対する事前分布や、データが従う確率分布を指定します。
stan
model {
  mu ~ normal(0, 1); // muの事前分布
  y ~ normal(mu, 1); // 観測データyの分布
}
  1. 生成量ブロック (generated quantities block)
    モデルから生成される追加の量を計算するためのブロックです。ここでは、例えば予測値や、モデルのフィットを評価するための統計量を計算します。このブロックで生成された値は、MCMCサンプルとともに出力されます。
stan
generated quantities {
  real y_pred = normal_rng(mu, 1); // 新しいデータポイントの予測
}

PythonからStanモデルを定義する方法

PythonからStanモデルを実行するには、pystanやcmdstanpyを使ってモデルを定義し、コンパイルする必要があります。ここでは、pystanを使用した方法を紹介します。

  1. Stanコードの定義
    Stanモデルは文字列としてPython内で定義されます。この文字列には上記で説明したStanのコードが含まれます。
Python
model_code = """
data {
  int<lower=0> N;
  real y[N];
}
parameters {
  real mu;
}
model {
  mu ~ normal(0, 1);
  y ~ normal(mu, 1);
}
generated quantities {
  real y_pred = normal_rng(mu, 1);
}
"""
  1. モデルのコンパイル
    次に、このStanコードをコンパイルします。コンパイルされたモデルはPythonオブジェクトとして保存され、後でデータを使ってサンプリングを行います。
Python
import pystan
stan_model = pystan.StanModel(model_code=model_code)
  1. データの準備
    コンパイルしたモデルにデータを入力するため、適切な形式でデータを準備します。データは辞書形式で渡されます。
Python
data = {'N': 10, 'y': [2.9, 3.1, 2.5, 2.8, 3.0, 3.2, 2.7, 3.0, 3.1, 2.9]}
  1. モデルの実行
    最後に、データをモデルに渡してサンプリングを行います。結果はfitオブジェクトに保存され、解析や可視化に使用されます。
Python
fit = stan_model.sampling(data=data, iter=1000, chains=4)

このようにして、PythonからStanモデルを定義し、データに基づいた推定を行うことができます。これで、基本的なStanモデルをPythonで実行する準備が整いました。

スポンサーリンク

Stanモデルの実行と結果の解析

Stanモデルを定義し、データを準備したら、次はモデルの実行と結果の解析を行います。この章では、Stanモデルのコンパイルとサンプリングのプロセス、結果の解析方法について詳しく解説します。

モデルのコンパイルとサンプリング

Stanを使用してベイズ推定を行うためには、まずStanモデルをコンパイルし、次にサンプリングを行う必要があります。サンプリングとは、StanがMCMC（マルコフ連鎖モンテカルロ）法を使ってパラメータの事後分布からサンプルを生成するプロセスです。

  1. Stanモデルのコンパイル
    まず、定義したStanコードをコンパイルします。コンパイルは一度行えば、以後はそのモデルを再利用できます。
Python
stan_model = pystan.StanModel(model_code=model_code)
  1. サンプリングの実行
    次に、データを使ってサンプリングを実行します。samplingメソッドを呼び出すと、指定したイテレーション数とチェーン数に基づいて、StanがMCMCサンプルを生成します。
Python
fit = stan_model.sampling(data=data, iter=1000, chains=4)
  • iterはMCMCの反復回数を指定し、chainsは独立して実行するチェーンの数を指定します。複数のチェーンを使用することで、サンプルが安定し、収束をチェックしやすくなります。

サンプルデータの準備とモデルへの投入

サンプルデータは、Stanモデルに入力するデータセットです。データはPythonの辞書形式で渡され、辞書のキーはStanのデータブロックで定義した変数名と一致させます。

例えば、以下のようにデータを準備し、Stanモデルに渡します。

Python
data = {
    'N': 10,
    'y': [2.9, 3.1, 2.5, 2.8, 3.0, 3.2, 2.7, 3.0, 3.1, 2.9]
}
fit = stan_model.sampling(data=data, iter=1000, chains=4)

Stanはこのデータを使って、指定されたモデルに基づき、muの事後分布をサンプリングします。

結果の解析と可視化

サンプリングが完了すると、Stanは推定されたパラメータの事後分布や、関連する統計量を返します。これをfitオブジェクトから取得し、解析や可視化を行います。

  1. 結果の表示
    fitオブジェクトのprintメソッドを使うと、推定されたパラメータの要約統計量が表示されます。ここには、パラメータの平均、標準誤差、95%信頼区間などが含まれます。
Python
print(fit)
  1. トレースプロットの作成
    トレースプロットは、MCMCサンプルの収束を確認するための重要なツールです。トレースプロットでは、各チェーンのサンプルがどのように変動しているかを視覚化します。matplotlibを使って簡単にプロットできます。
Python
import matplotlib.pyplot as plt
fit.plot()
plt.show()
  1. 事後分布の分析
    パラメータの事後分布を詳細に分析することで、モデルの結果を深く理解できます。事後分布は、推定されたパラメータがどの範囲に存在するかを示します。
Python
posterior_samples = fit.extract()
mu_samples = posterior_samples['mu']

plt.hist(mu_samples, bins=30, density=True)
plt.title('Posterior distribution of mu')
plt.xlabel('mu')
plt.ylabel('Density')
plt.show()

上記のコードでは、muの事後分布をヒストグラムとして表示します。これにより、muの推定値がどのように分布しているかを視覚的に確認できます。

  1. 他の解析手法
    Stanは、さらに高度な結果解析を行うためのさまざまな機能を提供しています。例えば、予測チェック（Posterior Predictive Check）を行ってモデルの適合性を評価したり、相関行列を分析してパラメータ間の依存関係を調べたりすることができます。
スポンサーリンク

Stanの高度な機能

Stanは基本的なベイズモデリングだけでなく、複雑な統計モデルを構築するための多くの高度な機能を提供しています。この章では、階層モデル、カスタム分布の定義、パフォーマンスの最適化など、Stanの高度な機能について詳しく説明します。

階層モデルと複雑なモデル構造

階層モデル（階層ベイズモデル）は、複数のレベルでデータが階層的に構造化されている場合に使用される統計モデルです。例えば、学校ごとに異なる生徒の成績を分析する際、学校レベルと生徒レベルの効果を別々にモデル化できます。

階層モデルでは、データに対して異なるレベルのパラメータを定義し、それらがどのように関連しているかをモデル化します。Stanでは、このような階層モデルを直感的に定義できます。

以下に、Stanでの階層モデルの例を示します。

stan
data {
  int<lower=0> J;  // グループの数（例えば、学校の数）
  vector[J] y;     // 各グループの観測値
  real<lower=0> sigma_y; // 観測誤差の標準偏差
}
parameters {
  real mu;          // 全体平均
  real<lower=0> tau; // グループ間の標準偏差
  vector[J] theta;   // 各グループの効果
}
model {
  theta ~ normal(mu, tau); // グループ効果の事前分布
  y ~ normal(theta, sigma_y); // 観測値の分布
}

このモデルでは、各グループ（学校）の効果thetaが全体平均muからどの程度離れているかをモデル化しています。tauはグループ間のばらつきを示し、sigma_yは観測誤差を表します。

カスタム分布の定義

Stanには多くの標準的な確率分布が組み込まれていますが、特殊な分布や問題に特化した分布が必要な場合、自分でカスタム分布を定義することもできます。カスタム分布の定義には、Stanのfunctionsブロックを使用します。

例えば、ゼロインフレーションポアソン分布をカスタム分布として定義する場合のStanコードは次のようになります。

stan
functions {
  real zero_inflated_poisson_lpmf(int y, real lambda, real zi) {
    if (y == 0)
      return log_sum_exp(bernoulli_lpmf(1 | zi), bernoulli_lpmf(0 | zi) + poisson_lpmf(0 | lambda));
    else
      return bernoulli_lpmf(0 | zi) + poisson_lpmf(y | lambda);
  }
}
data {
  int<lower=0> N;
  int<lower=0> y[N];
}
parameters {
  real<lower=0> lambda;
  real<lower=0, upper=1> zi;
}
model {
  for (n in 1:N)
    y[n] ~ zero_inflated_poisson(lambda, zi);
}

このカスタム分布は、ゼロが過剰に観測されるデータをモデル化するために使用されます。

パフォーマンスの最適化とヒント

Stanのパフォーマンスは、モデルの複雑さやデータのサイズに大きく依存します。モデルの実行速度を最適化するために、いくつかのヒントがあります。

  1. 初期値の指定
    初期値を適切に指定することで、MCMCサンプルがより速く収束することがあります。特に、パラメータのスケールが大きく異なる場合や、事前に信頼できる初期値がわかっている場合には効果的です。
Python
init_values = {'mu': 0.5, 'tau': 1.0}
fit = stan_model.sampling(data=data, iter=1000, chains=4, init=[init_values]*4)
  1. 再パラメータ化
    モデルの再パラメータ化（リファクタリング）は、サンプリングの効率を大幅に向上させることがあります。例えば、非中心化されたパラメータ化（non-centered parameterization）は、特に階層モデルで効果的です。

非中心化されたパラメータ化の例:

stan
parameters {
  real mu;
  real<lower=0> tau;
  vector[J] theta_raw;
}
transformed parameters {
  vector[J] theta = mu + tau * theta_raw;
}
model {
  theta_raw ~ normal(0, 1);
  y ~ normal(theta, sigma_y);
}
  1. 並列化
    cmdstanpyを使用することで、Stanモデルのサンプリングを複数のCPUコアで並列化することができます。これにより、複数のチェーンを並行して実行し、サンプリングの時間を短縮できます。
Python
fit = stan_model.sample(data=data, chains=4, parallel_chains=4)
  1. プリコンパイル済みモデルの再利用
    モデルを頻繁に実行する場合、プリコンパイル済みのStanモデルを再利用することで、コンパイル時間を節約できます。cmdstanpyでは、コンパイル済みのモデルを保存し、再利用することが容易です。
スポンサーリンク

実例: PythonでStanを使った実践モデリング

ここでは、具体的なケーススタディを通じて、PythonでStanを使ったベイズモデリングの実際のプロセスを紹介します。例として、ベイズ回帰モデルを構築し、データを使って推定を行います。

ケーススタディ: ベイズ回帰モデル

回帰分析は、ある変数（目的変数）が他の変数（説明変数）とどのように関係しているかをモデル化する手法です。ベイズ回帰では、この関係性をベイズ推定のフレームワーク内で表現し、パラメータの不確実性を考慮に入れたモデルを構築します。

今回は、単純な線形回帰モデルを例にとり、Stanを使ってそのパラメータを推定してみましょう。

データの準備

まず、回帰モデルに使用するデータを準備します。ここでは、簡単な人工データを生成します。

Python
import numpy as np
import matplotlib.pyplot as plt

# データ生成
np.random.seed(42)
N = 100  # サンプル数
x = np.random.normal(0, 1, N)
y = 2.5 * x + np.random.normal(0, 1, N)

# データをプロット
plt.scatter(x, y)
plt.xlabel('x')
plt.ylabel('y')
plt.title('Synthetic Data for Bayesian Linear Regression')
plt.show()

このデータは、yがxに対して線形関係を持つものです。回帰係数は2.5ですが、ノイズが加わっています。

Stanモデルの定義

次に、このデータにフィットするStanモデルを定義します。線形回帰モデルでは、目的変数yは説明変数xに対して線形な関係を持ち、パラメータとして回帰係数と切片を持ちます。

Python
model_code = """
data {
  int<lower=0> N;         // サンプル数
  vector[N] x;            // 説明変数
  vector[N] y;            // 目的変数
}
parameters {
  real alpha;             // 切片
  real beta;              // 回帰係数
  real<lower=0> sigma;    // 残差の標準偏差
}
model {
  y ~ normal(alpha + beta * x, sigma);  // 線形モデル
}
"""

このStanモデルでは、yはxに対してalpha + beta * xという線形関係を持ち、残差（誤差）は正規分布に従うと仮定しています。Stanでは、alpha（切片）とbeta（回帰係数）、およびsigma（誤差の標準偏差）を推定します。

モデルのコンパイルとサンプリング

次に、定義したStanモデルをコンパイルし、データを使ってサンプリングを行います。

Python
import pystan

# データの準備
data = {'N': N, 'x': x, 'y': y}

# モデルのコンパイル
stan_model = pystan.StanModel(model_code=model_code)

# サンプリングの実行
fit = stan_model.sampling(data=data, iter=1000, chains=4)

# 結果の表示
print(fit)

このコードを実行すると、alpha、beta、sigmaの推定値とその不確実性（事後分布）が表示されます。

結果の解釈と可視化

サンプリングの結果を解釈し、可視化することで、モデルがデータにどのように適合しているかを確認します。

Python
# パラメータの事後分布を取得
posterior_samples = fit.extract()
alpha_samples = posterior_samples['alpha']
beta_samples = posterior_samples['beta']

# 回帰線のプロット
plt.scatter(x, y, color='blue', alpha=0.5)
plt.plot(x, alpha_samples.mean() + beta_samples.mean() * x, color='red')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Fitted Line using Bayesian Linear Regression')
plt.show()

このプロットでは、青い点が観測データ、赤い線が推定された回帰直線を示しています。これにより、モデルがデータにうまくフィットしているかを視覚的に確認できます。

また、パラメータの事後分布をヒストグラムとして表示し、それぞれのパラメータがどの範囲にあるかを視覚化することも重要です。

Python
plt.hist(alpha_samples, bins=30, density=True)
plt.title('Posterior distribution of alpha')
plt.xlabel('alpha')
plt.ylabel('Density')
plt.show()

plt.hist(beta_samples, bins=30, density=True)
plt.title('Posterior distribution of beta')
plt.xlabel('beta')
plt.ylabel('Density')
plt.show()

これらのヒストグラムにより、alphaやbetaの不確実性や信頼区間を確認できます。

このケーススタディでは、Stanを使って単純なベイズ回帰モデルをPythonで実装し、データに対してモデルをフィットさせました。このようにして、Stanを利用することで、複雑な統計モデルを効率的に構築し、その結果を直感的に解釈することが可能です。

スポンサーリンク

まとめ

PythonでStanを使用することで、複雑なベイズ統計モデルを効率的に構築し、データに対する推論や予測を行うことができます。本記事では、Stanの基本的な使い方から高度な機能までをカバーし、実際のケーススタディを通じてStanの実用性を示しました。

PythonでStanを使用するメリットと限界

Pythonを介してStanを利用することには多くのメリットがあります。Pythonの広範なライブラリとの連携により、データの前処理、モデルの実行、結果の解析、そして可視化までを一貫して行えるのは大きな強みです。特に、pystanやcmdstanpyといったライブラリを使えば、Stanのパワフルなモデリング機能を手軽に活用できます。

一方で、Stanの使用にはいくつかの限界も存在します。特に複雑なモデルや大規模なデータセットを扱う場合、サンプリングに要する時間が長くなることがあります。また、MCMC法の特性上、収束の確認や適切な初期値の設定が必要です。

今後の展望と参考資料

Stanは継続的にアップデートされており、新しい機能や最適化手法が追加されています。特に、最新のバージョンでは、さらなる性能向上や新しい分布、モデル構造のサポートが追加されており、より広範なモデリングが可能になっています。

また、コミュニティも活発で、多くのチュートリアルや参考資料がオンラインで提供されています。これらの資料を活用することで、Stanの知識をさらに深め、より高度なモデリング技術を習得することができます。

参考として、以下の公式ドキュメントやコミュニティリソースを参照することをお勧めします。

brian
brian

ここまで読んでいただきありがとうございます！

UdemyのPythonコースにはオンラインで学習ができる動画コンテンツがたくさんあります。

当ブログのような文章メインの説明では足りない箇所を補えると思うので、もっと詳しく勉強したいという方はぜひチェックしてみてください！

コメント

タイトルとURLをコピーしました