【Python】 PystanによるMCMCを用いた回帰分析

はじめに

　回帰分析(フィッティング)の方法として、最小二乗法やカイ二乗検定、最尤法が一般的な方法としてあり、それらは損失関数を最小化/最大化することで、データに合うモデルパラメーターを決めている。ただ、パラメータが増えたり関数が複雑になると、最適な解を求めるには数値計算が難しくなり、近似的に解を求める必要がある。そのような近似解を求める方法の一つとして、今回はマルコフ連鎖モンテカルロ法を紹介する。

• マルコフ連鎖モンテカルロ法(MCMC:Markov chain Monte Carlo methods)
  事後分布を数値計算から求めることが困難な場合、事前分布と事後分布にマルコフ連鎖を仮定し、モンテカルロ法とベイズの定理により確率的に事後分布(次の状態)を計算し、近似値を探索する方法。

• 変分推論(VI:Variational Inference)
  近似分布と事後分布のロス関数をKLダイバージェンスで定義し、ロス関数を勾配法などを用いて最小化することで近似値を探索する方法。

ChatGPTに聞いてみた

MCMC (マルコフ連鎖モンテカルロ) は、確率分布からサンプリングするアルゴリズムのファミリーです。これらのアルゴリズムは、確率の推定、最適化問題の解決、およびベイジアン推論の実行に使用されます。

MCMC の背後にある基本的な考え方は、定常分布として目的の確率分布を持つマルコフ連鎖を構築することです。次に、アルゴリズムは、サンプルが定常分布に収束することを期待して、マルコフ連鎖から一連のランダム サンプルを生成します。このアルゴリズムは、提案分布を使用して候補サンプルを生成し、現在のサンプルに対する尤度に基づいて候補サンプルを受け入れるか拒否します。時間の経過とともに、マルコフ連鎖はサンプル空間のさまざまな領域を訪れ、連鎖によって生成されたサンプルは目的の分布に収束します。

MCMC アルゴリズムは、機械学習、計算生物学、経済学など、さまざまな分野で広く使用されています。 MCMC アルゴリズムの一般的な例には、Metropolis-Hastings、Gibbs サンプリング、ハミルトニアン モンテカルロなどがあります。

MCMC アルゴリズムには、高次元分布を処理する機能やベイジアン推論を実行する機能など、いくつかの利点があります。ただし、MCMC アルゴリズムは収束が遅くなる可能性があり、ミキシングが不十分である、自己相関が高い、提案分布の選択に対する感度などの問題が発生する可能性があります。これらの課題を克服するには、当面の問題に適した MCMC アルゴリズムを選択し、アルゴリズムのパラメーターを慎重に調整して良好なパフォーマンスを確保することが重要です。

フィッティングにおけるベイズ統計の解釈

　あるデータ$x$のとき、$y(x,\theta)$の分布が観測値として求まる。その分布は自然に正規分布になるとする。そのとき、観測値の分布を最も説明する分布$P(x|\theta)$は$\theta$のみに依存する。$\theta$は正規分布の平均値と分散$\mu, \sigma$であり、$y(x|\theta)$である。つまり観測値に最も合うパラメータ$\theta$を求めることができる。

Pystan

PystanはStanのPythonライブラリである。Stanは統計モデリングのための最先端のオープンソースプラットフォームで、C++で提供されている。

Pythonコード

import os
import sys
import numpy as np
import matplotlib.pyplot as plt
from multiprocessing import Process, freeze_support
import pymc as pm
import arviz as az
import csv

RANDOM_SEED = 8927
rng = np.random.default_rng(RANDOM_SEED)

# サンプルデータ生成
N = 100
alpha_real = 2.5
beta_real = 0.9
gamma_real = 3.0
eps_real = np.random.normal(0, 50, N)

# モデル: y = α + βx + γxx + Noise
x = np.random.normal(10, 3, N)
x = np.sort(x)
y = alpha_real + beta_real*x + gamma_real*x*x + eps_real

if __name__ == '__main__':
    with pm.Model() as model:
        # 1. 事前分布(提案分布) : 観測されていない確率分布
        alpha = pm.Normal('alpha', mu=100, sigma=10) # αには平均0, 標準偏差10の正規分布を仮定
        beta = pm.Normal('beta', mu=100, sigma=10) # βには平均0, 標準偏差1の正規分布を仮定
        gamma = pm.Normal('gamma', mu=100, sigma=10) # γには平均0, 標準偏差1の正規分布を仮定
        epsilon = pm.HalfCauchy('sd', beta=5.0) # 標準偏差の分布には最頻値5の半コーシー分布を仮定、裾野が広いため、たまに大きな値が出る

        # 2. 尤度の計算 : 観測された確率分布、引数にobservedが存在すれば、尤度関数として定義される
        mu = alpha + beta*x + gamma*x*x
        y_pred = pm.Normal('y_pred', mu=mu, sigma=epsilon, observed=y) # 
        print('\nモデルの基本パラメーター\n',model.basic_RVs, '\nフリーパラメーター\n', model.free_RVs, '\n観測パラメーター\n', model.observed_RVs,'\n')

        # 3. MCMCの実行
        step = pm.NUTS() #MCMC モデル BinaryGibbsMetropolis', 'BinaryMetropolis', 'CategoricalGibbsMetropolis', 'CauchyProposal', 'CompoundStep', 'DEMetropolis', 'DEMetropolisZ', 'HamiltonianMC', 'LaplaceProposal', 'Metropolis', 'MultivariateNormalProposal', 'NUTS', 'NormalProposal', 'PoissonProposal', 'STEP_METHODS', 'Slice', 'UniformProposal'] 
        trace = pm.sample(draws=3000, tune=1000, step=step, chains=2, cores=32, return_inferencedata=True)
        trace.to_netcdf("data/trace.nc")
        print('\nCorrectly Done : MCMC Sampling!')

        # 4. 結果の確認
        print('\n', pm.summary(trace))
        pm.plot_trace(trace)
        fig = plt.gcf()
        fig.savefig('mcmc_sampling.jpg')
        print('\nCorrectly Done : Save Sampling Result PLot!')
        ppc = pm.sample_posterior_predictive(trace, model=model, random_seed=rng, return_inferencedata=False)
        alpha_m = trace['posterior']['alpha'].mean().values
        beta_m = trace['posterior']['beta'].mean().values
        gamma_m = trace['posterior']['gamma'].mean().values
        print('\nalpha mean', alpha_m,'\nbeta mean', beta_m,'\ngamma mean', gamma_m,'\n')
        fig = plt.figure()
        plt.plot(x, y, '.')
        plt.plot(x, alpha_m + beta_m*x + gamma_m*x*x)
        idx = np.argsort(x)
        x_ord = x[idx]
        sig0 = az.hdi(ppc['y_pred'], hdi_prob=0.5)[idx]
        sig1 = az.hdi(ppc['y_pred'], hdi_prob=0.95)[idx]
        plt.fill_between(x_ord, sig0[:,0], sig0[:,1], color='gray', alpha=1)
        plt.fill_between(x_ord, sig1[:,0], sig1[:,1], color='gray', alpha=0.5)
        fig.savefig('mcmc_fit.jpg')
        print('\nCorrectly Done : Save fitting result PLot!')

実行結果

平均値がα=68.2、β=-14.1、γ=3.8と結構はなれた数値になったが、フィッティング画像を見るによくフィットしているように見える。

まとめ

今回は、PystanによるMCMCを紹介した。