気候予測データの解析環境を構築する。その7―MCMCサンプラーの実行環境を構築する。

はじめに

気候予測データは、将来の平均的な「気候」状態を予測するためのもので、ある地域の気温・降水量の平均値や変動の標準偏差などの統計量がどう変化するか？を予測しています。(ココが知りたい地球温暖化より)

www.cger.nies.go.jp

気温など、連続した数値でかつ正規分布に従うことがわかっていれば、平均値や分散を求めることは簡単です。しかし、ガンマ分布に従う降水量や離散値である猛暑日日数の変化を単純な算術平均で表すことには問題があります。そのようなデータを扱うには"正規分布以外の確率分布を扱える統計モデル"が有用です。

“階層ベイズモデル”は、様々な確率分布を扱い複雑な統計モデルを構築することができ、マルコフ連鎖モンテカルロ法を用いてサンプリングをすることでパラメータを推計することができます。

Pythonでマルコフ連鎖モンテカルロ法を使うためのライブラリ（MCMCサンプラー）には、PyMC、Pyro, NumPyro, TensorFlowProbability, PyStan, PyJAGSなどがありますが、本研究室では"PyMC"を使っています。

www.pymc.io

ということで、PyMCの実行環境を構築していきます。

PyMCのインストール

Python標準の"venv"を使って、MCMC用の環境を作成し、PyMCをインストールします。

# 最新の状態に更新　
sudo apt update
sudo apt upgrade

# venvのインストール
sudo apt install python3.10-venv

# 仮想環境の構築
python3 -m venv v-MCMC

# v-MCMCに入る →プロンプトの頭に（v-MCMC）と表示される
source v-MCMC/bin/activate

# ライブラリの追加
pip3 install scipy pandas
pip3 install plotly kaleido matplotlib
pip3 install ipykernel
pip3 install  cloudpickle cachetools arviz
pip3 install fastprogress typing-extensions
pip3 install pytensor

# PyMCのインストール
pip3 install pymc
pip3 install git+https://github.com/pymc-devs/pymc-experimental.git

VScodeからv-MCMC環境を呼び出す

↓を参考に、VSCodeにvenv環境を認識してもらいます。

qiita.com

GUIでは、右下の「管理」から「設定」を呼び出し、

"venv"で検索して「項目の追加」をクリック

"v-MCMC"と入力してOK

追加される

「表示」→「コマンドパレット」からJupyterインタラクティブシェルを起動

右上のカーネル名をクリック→「別のカーネルを選択」

「Python環境」を選択

「v-MCMC」を選択

カーネルが変更される

一般化線形モデルを動かしてみる

PyMCの公式には様々なチュートリアルが公開されています。動作確認に“GLM: Linear regression”を実行してみます。

https://www.pymc.io/projects/docs/en/latest/learn/core_notebooks/GLM_linear.html

まずは環境設定と擬似データを作成

# ライブラリの読み込み
import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymc as pm
import xarray as xr
from pymc import HalfCauchy, Model, Normal, sample

# 環境設定
RANDOM_SEED = 8927
rng = np.random.default_rng(RANDOM_SEED)

%config InlineBackend.figure_format = 'retina'
az.style.use("arviz-darkgrid")

# 擬似データ作成
size = 200
true_intercept = 1
true_slope = 2

x = np.linspace(0, 1, size)
# y = a + b*x の線形回帰モデル
true_regression_line = true_intercept + true_slope * x
# データにノイズを加える
y = true_regression_line + rng.normal(scale=0.5, size=size)

data = pd.DataFrame(dict(x=x, y=y))

# データのグラフを描いてみる
fig = plt.figure(figsize=(7, 7))
ax = fig.add_subplot(111, xlabel="x", ylabel="y", title="Generated data and underlying model")
ax.plot(x, y, "x", label="sampled data")
ax.plot(x, true_regression_line, label="true regression line", lw=2.0)
plt.legend(loc=0);

次にモデルを組み立てます。チュートリアルでは一気にサンプリングまでしていますが、ここではまず、モデルを組み立てるところまで。

with Model() as model:  # model specifications in PyMC are wrapped in a with-statement
    # Define priors
    sigma = HalfCauchy("sigma", beta=10)
    intercept = Normal("Intercept", 0, sigma=20)
    slope = Normal("slope", 0, sigma=20)

    # Define likelihood
    likelihood = Normal("y", mu=intercept + slope * x, sigma=sigma, observed=y)

チュートリアルにはありませんが、以下のコマンドでモデルの内容を表示できます。

# 定義した各パラメータの情報
display(model.model)

# モデル構造
modeldag = pm.model_to_graphviz(model)
display(modeldag)

サンプリングを実行

with model:
    idata = sample(3000)

結果を表示します。 ※チュートリアルでは"bambi"などなるものを使って、Rで使うようなチルダ式を使ってモデルを作成する方法が紹介されていますが、ここでは割愛します。

# サンプリング結果
az.plot_trace(idata, figsize=(10, 7));

# 予測結果を計算
idata.posterior["y_model"] = idata.posterior["Intercept"] + idata.posterior["slope"] * xr.DataArray(x)

# 予測結果を描く
_, ax = plt.subplots(figsize=(7, 7))
az.plot_lm(idata=idata, y="y", num_samples=100, axes=ax, y_model="y_model")
ax.set_title("Posterior predictive regression lines")
ax.set_xlabel("x");

動きましたか？

次回はこのPyMCを使って、 「日本域CMIP6データ」を解析したいと思います。