見出し画像

時系列分析入門【第４章 その１】状態空間モデル：ローカルレベルモデルをPythonとPyMC Ver.5 で実践する

ネイピア DS

この記事は、テキスト「RとStanではじめる 心理学のための時系列分析入門」の第４章「RStanによる状態空間モデル」のRスクリプト・Stanコードをお借りして、PythonとPyMC Ver.5 で「実験的」に実装する様子を描いた統計ドキュメンタリーです。

取り扱いテーマは「ローカルレベルモデル」です。
状態空間モデルの基礎編に相当します。
３つの方法でベイズ推論を実施します。

  • pymcのGaussianRandomWalk

  • pymcのAR

  • pymc_experimental.statespace

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

テキストで使用するデータは、R・Stan等のサンプルスクリプトで指定されています。
サンプルスクリプトは著者GitHubサイトからダウンロードして取得できます。



第４章 RStanによる状態空間モデル 

この記事は「4.1 RStanによる状態空間モデル」を実践いたします。

インポート（＋環境準備）

この記事で用いるライブラリをインポートします。

### インポート

# 基本
import numpy as np
import pandas as pd

# PyMC
import pymc as pm
from pymc_experimental.statespace import structural as st
import arviz as az

# グラフ描画
import matplotlib.pyplot as plt
plt.rcParams['font.family'] = 'Meiryo'

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

【補足】PyMC環境の構築
AnacondaでPyMC環境を構築する方法を次のリンクの「pymc環境の構築」の箇所に記載しています。
ご参考まで。

【補足】pymc_experimentalのインストール
pymc_experimentalの公式サイトに掲載のインストール方法は、Gitリポジトリからインストールする方法です。
そこで、まずGitアプリをインストールして、次にpymc_experimentalをインストールします。

■ Gitアプリのインストール
私はこちらのサイトからGitアプリをダウンロードしました。

Windows版のインストールの手順は、こちらのサイトの情報を参考にしました。

■pymc_experimentalのインストール
pymc_experimental公式サイトのインストールコードをそのまま実行しました。

pip install git+https://github.com/pymc-devs/pymc-experimental.git


ローカルレベルモデルの概要

ローカルレベルモデルの数式表現です。

$${\mu_t = \mu_{t-1} + u_t, \quad u_t \sim \text{Normal}\ (0,\ \sigma^2_{u})}$$
$${y_t = \mu_t + v_t, \quad v_t \sim \text{Normal}\ (0,\ \sigma^2_{v})}$$

テキストより引用

１つ目の数式を「状態方程式」「システムモデル」と呼び、２つ目の式を「観測方程式」「観測モデル」と呼びます。
$${y_t}$$は時間$${t}$$の観測値です。観測ノイズ$${v_t}$$を含みます。
$${\mu_t}$$は時間$${t}$$で観測できない「状態」を示す値です。１期前の時間$${t-1}$$の状態$${\mu_{t-1}}$$にノイズ$${u_t}$$を加味しています。

ローカルレベルモデルの観測値と状態の生成過程は正規分布と仮定されており、次のように表現できます。

$${\mu_t \sim \text{Normal}\ (\mu_{t-1},\ \sigma^2_u)}$$
$${y_t \sim \text{Normal}\ (\mu_t,\ \sigma^2_v)}$$

テキストより引用

また、ローカルレベルモデルの状態は「ランダムウォーク」です。

【ランダムウォーク】
$${y_t = y_{t-1} + \varepsilon_t, \quad \varepsilon_t \sim \text{Normal}\ (0,\ \sigma^2)}$$

テキストより引用

図示します。

■ 小まとめ
ランダムウォークと正規分布に従う確率変数であることを手がかりにして、最初に、pymcのガウスランダムウォークおよび自己回帰モデルを検討します。
その後、pymc_experimentalの状態空間モデルを検討します。


データの準備

２つのデータを読み込みます。
テキストのRコードに記述されたデータをお借りします。

① 基本モデル用データ：欠損値の無いデータ
テキストの「4.1.4 RStanの使い方」で用いられるデータです。
Rコードのデータをお借りして、pandasのデータフレーム df に格納します。

### データ作成１

# 体重仮想データ
Y = [69.9, 69.3, 69.8, 70.4, 70.6, 70.3, 69.9, 69.6, 69.1, 69.5,
     70.1, 70.3, 71.0, 69.9, 69.8, 70.1, 69.8, 69.7, 69.3, 69.1]

# データフレーム化
df = pd.DataFrame({'Day': range(1, len(Y)+1), 'Weight': Y})

# 描画 ※図4.3に相当
plt.figure(figsize=(10, 4))
plt.plot(df['Day'], df['Weight'])
plt.xlabel('Day')
plt.ylabel('Weight(kg)')
plt.xticks([5, 10, 15, 20])
plt.grid(lw=0.2)
plt.show()

【実行結果】
20日間の体重は69kg～71kgを上下しています。

② 将来予測モデル用データ：欠損値のあるデータ
テキストの「4.1.5 欠損値の補間と将来の予測は同じ」で用いられています。
データフレーム df2 に格納します。
あわせて１０期先予測用のデータ y_fore, idx_fore を作成します。

### データ作成２

# 体重仮想データ
Y2 = [69.9, 69.3, 69.8, 70.4, 70.6, 70.3, np.nan, 69.6, 69.1, 69.5,
      70.1, 70.3, 71.0, np.nan, 69.8, 70.1, 69.8, 69.7, 69.3, 69.1]

# データフレーム化
df2 = pd.DataFrame({'Day': range(1, len(Y2)+1), 'Weight': Y2})
display(df2)

# PyMCモデル用データの作成（予測期間を考慮）
# 予測期間
pred_period = 10

# 観測データ＋予測期間データの生成
y_fore = np.concatenate([df2['Weight'].values, np.repeat(np.nan, pred_period)])
idx_fore = np.concatenate([df2.index, np.arange(df2.index[-1]+1,
                                                df2.index[-1]+pred_period+1)])

【実行結果】
７日目、１４日目に欠損値が含まれます。

将来予測モデル用データ df2 を描画します。

### 描画
plt.figure(figsize=(10, 4))
plt.plot(df2['Day'], df2['Weight'])
plt.xlabel('Day')
plt.ylabel('Weight(kg)')
plt.xticks([5, 10, 15, 20])
plt.grid(lw=0.2)
plt.show()

【実行結果】
折れ線グラフが途切れています！


ガウスランダムウォーク

状態方程式にpymcのGaussianRandomWalkを用いるモデルです。
次のリンクでPyMC公式サイトのGaussianRandomWalkのページを開けます。

１．基本モデル

① モデルの定義
状態「mu」 に pm.GaussianRandomWalk()を用いています。

### モデルの定義

# モデルの定義
with pm.Model() as model1_1:
    
    ### データ関連定義
    # coordの定義
    model1_1.add_coord('data', values=df.index, mutable=True)
    # dataの定義
    y = pm.MutableData('y', value=df['Weight'].values, dims='data')
    
    ### 事前分布
    sigmaU = pm.Uniform('sigmaU', lower=0, upper=100)
    sigmaV = pm.Uniform('sigmaV', lower=0, upper=100)
    
    ### muの計算
    init_dist = pm.Normal.dist(mu=0, sigma=100)
    mu = pm.GaussianRandomWalk('mu', sigma=sigmaU, init_dist=init_dist,
                               dims='data')
    
    ### 尤度
    likelihood = pm.Normal('likelihood', mu=mu, sigma=sigmaV, observed=y,
                           dims='data')

【実行結果】出力なし

② モデルの表示・可視化

### モデルの表示
model1_1

【実行結果】
RandomWalkの引数の１つ目は指定した初期分布、２つ目は通常？の分布であり、平均０、標準偏差sigmaUの正規分布です。

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

【実行結果】
シンプルなモデルです。

③ 事後分布からのサンプリング
MCMCを実行します。
NUTSサンプラーにnumpyroを使用するので爆速です。
実行時間は約12秒でした。

### 事後分布からのサンプリング ※NUTSサンプラーにnumpyroを使用　12秒
with model1_1:
    idata1_1 = pm.sample(draws=2000, tune=2000, chains=4, target_accept=0.9,
                         nuts_sampler='numpyro', random_seed=1234)

【実行結果（一部分）】

④ サンプリングデータの確認
$${\hat{R}}$$、事後分布の要約統計量、トレースプロットを確認します。
事後分布の収束確認はテキストにならって$${\hat{R} \leq 1.1}$$としています。

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

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

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

パラメータ等の事後分布の要約統計量です。
テキストの表4.1に相当します。

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

【実行結果】
状態「mu」の時間インデックス 0～19（時間 1～20 に相当）、および、状態ノイズsigmaU、観測ノイズsigmaVの要約統計量を表示しています。
テキストの表4.1の平均：mean、標準偏差sdの値に近いのではないでしょうか。
なお、2.5%・97.5%に関して、この図は95%HDI区間を示しています。
一方でテキストは確信区間（信用区間）を示しており、計算基礎が異なっています。

トレースプロットを確認しましょう。

### トレースプロットの表示
pm.plot_trace(idata1_1, combined=True, figsize=(12, 8))
plt.tight_layout();

【実行結果】
収束している感じがします。
若干、発散のバーコードが見られます。

⑤ サンプリングデータを描画
状態「mu」の事後平均値と95%HDIを描画します。
テキストの図4.4に相当します。
同じ処理を何度か実行しますので、描画処理を関数化いたします。

### muの事後分布の平均値、50%・95%HDI区間の描画関数の定義

def plot_mu_posterior(idata):
    # 事後平均・95%HDIの算出
    tmp = idata.posterior
    mean_post = tmp.mu.stack(sample=('chain', 'draw')).mean(axis=1).data
    hdi_post95 = az.hdi(tmp, hdi_prob=0.95).mu.data
    hdi_post50 = az.hdi(tmp, hdi_prob=0.50).mu.data

    # 描画
    fig, ax = plt.subplots(figsize=(10, 4))
    # 観測値の描画
    ax.scatter(df['Day'], df['Weight'], color='blue', s=30, label='観測値')
    # 事後平均の描画
    ax.plot(df['Day'], mean_post, color='tab:blue', label='事後平均')
    # 95%HDI区間の描画
    ax.fill_between(df['Day'], hdi_post95[:, 0], hdi_post95[:, 1],
                    color='tab:blue', alpha=0.15, label='HDI 95%')
    # 50%HDI区間の描画
    ax.fill_between(df['Day'], hdi_post50[:, 0], hdi_post50[:, 1],
                    color='tab:blue', alpha=0.4, label='HDI 50%')
    # 修飾
    ax.set(xlabel='Day', ylabel='μの事後分布(kg)')
    ax.grid()
    ax.legend()
    plt.show()

描画処理を実行します。

### muの事後分布の平均値、50%・95%HDI区間の描画 ※図4.4に相当
plot_mu_posterior(idata1_1)

【実行結果】
テキストと近似しているように感じます（ホッとしました）。

続いて、plot_traceを用いて、sigmaUの事後分布と収束の様子を描画します。
テキストの図4.5と図4.6に相当します。

### sigma_Uの事後分布と収束の様子の描画 図4.5、4.6に相当
pm.plot_trace(idata1_1, var_names=['sigmaU'], combined=False, compact=False, 
              figsize=(10, 3), legend=True);

【実行結果】
バーンイン期間（tune）の描画は断念しましょう。

⑥ 事後予測
目的変数「Weight」の事後予測を行います。
テキストには無い処理です。

### 事後予測のサンプリングの実施
with model1_1:
    idata1_1.extend(pm.sample_posterior_predictive(idata1_1))

【処理結果】省略

目的変数「Weight」の事後予測サンプリングデータを描画します。
描画処理を関数化します。

### 目的変数の予測値の平均値、50%・95%HDI区間の描画関数の定義

def plot_post_predict(idata):
    # 事後平均・95%HDIの算出
    tmp = idata.posterior_predictive
    mean_post = tmp.likelihood.stack(sample=('chain', 'draw')).mean(axis=1).data
    hdi_post95 = az.hdi(tmp, hdi_prob=0.95).likelihood.data
    hdi_post50 = az.hdi(tmp, hdi_prob=0.50).likelihood.data

    # 描画
    fig, ax = plt.subplots(figsize=(10, 4))
    # 観測値の描画
    ax.scatter(df['Day'], df['Weight'], color='blue', s=30, label='観測値')
    # 事後平均の描画
    ax.plot(df['Day'], mean_post, color='tab:blue', label='予測平均')
    # 95%HDI区間の描画
    ax.fill_between(df['Day'], hdi_post95[:, 0], hdi_post95[:, 1],
                    color='tab:blue', alpha=0.15, label='HDI 95%')
    # 50%HDI区間の描画
    ax.fill_between(df['Day'], hdi_post50[:, 0], hdi_post50[:, 1],
                    color='tab:blue', alpha=0.4, label='HDI 50%')
    # 修飾
    ax.set(xlabel='Day', ylabel='Yの事後予測(kg)')
    ax.grid()
    ax.legend()
    plt.show()

描画処理を実行します。

### 目的変数の予測値の平均値、50%・95%HDI区間の描画
plot_post_predict(idata1_1)

【実行結果】
先の状態「mu」の描画と比べて、95%HDI区間は広く感じます。
観測ノイズをのせているからでしょうか。

２．将来予測モデル
欠損値の補間と１０期先予測を実施します。
テキストの図4.7、図4.9に相当する図を作成します。

将来予測モデル用に生成したデータ y_fore, idx_fore を用います。
これらのデータは、20日間の観測値（２日分は欠損値）に10日間の欠損値を追加したものです。
将来分も欠損値なのです。

① モデルの定義
欠損値を含むデータをモデル内の「data」で定義できません。
観測方程式（尤度）を定義する正規分布の引数 observed で観測データ y_fore を渡します。

### モデルの定義
# ※欠損値のデータはdataを利用できない。observedに目的変数を直接代入する

# モデルの定義
with pm.Model() as model1_2:
    
    ### データ関連定義
    # coordの定義
    model1_2.add_coord('data', values=idx_fore, mutable=True)
    # # dataの定義
    # y = pm.MutableData('y', value=y_4_pred, dims='data')
    
    ### 事前分布
    sigmaU = pm.Uniform('sigmaU', lower=0, upper=100)
    sigmaV = pm.Uniform('sigmaV', lower=0, upper=100)
    
    ### muの計算
    init_dist = pm.Normal.dist(mu=0, sigma=100)
    mu = pm.GaussianRandomWalk('mu', sigma=sigmaU, init_dist=init_dist,
                               dims='data')
    
    ### 尤度
    likelihood = pm.Normal('likelihood', mu=mu, sigma=sigmaV, observed=y_fore,
                           dims='data')

【実行結果】出力なし

② モデルの可視化

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

【実行結果】
12個の欠損値を推論する likelihood_unobserved が追加されました。
なお、モデルの表示「model1_2」を実行するとエラーになりました。

③ 事後分布からのサンプリング
MCMCを実行します。
NUTSサンプラーにnumpyroを指定するとエラーになりました。
欠損値を含むデータをnumpyroで利用する方法が分からないので、PyMC標準のサンプラーを用いることにします。
実行時間は約7分30秒でした。

### 事後分布からのサンプリング ※NUTSサンプラーにnumpyroを未使用　7分30秒
# 欠損値のある場合のnumpyroの使い方がわからない

with model1_2:
    idata1_2 = pm.sample(draws=2000, tune=2000, chains=4, target_accept=0.9,
                         random_seed=1234)

【実行結果】省略

④ サンプリングデータの確認
$${\hat{R}}$$、事後分布の要約統計量、トレースプロットを確認します。
事後分布の収束確認はテキストにならって$${\hat{R} \leq 1.1}$$としています。

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

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

【実行結果】
$${\hat{R}>1.1}$$のパラメータは「0」件です。
全てのパラメータが$${\hat{R} \leq 1.1}$$であることを確認できました。
また「likelihood_unobserved」および「likelihood」の推論値が追加されました。

パラメータ等の事後分布の要約統計量です。
テキストの表4.1に相当します。

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

【実行結果】
雰囲気をお楽しみ下さい。

トレースプロットを確認しましょう。

### トレースプロットの表示
pm.plot_trace(idata1_1, combined=True, figsize=(12, 8))
plt.tight_layout();

【実行結果】
収束している感じがします。
発散のバーコードが見られます。

⑤ サンプリングデータを描画
将来期間を含む状態「mu」の事後平均値と95%HDIを描画します。
テキストの図4.9に相当します。
同じ処理を何度か実行しますので、描画処理（予測用）を関数化いたします。

### muの将来予測値の平均値、50%・95%HDI区間の描画関数の定義
def plot_mu_posterior2(idata):
    
    ## 事後平均・95%HDI・50%HDIの算出
    tmp = idata.posterior
    mean_post = tmp.mu.stack(sample=('chain', 'draw')).mean(axis=1).data
    hdi_post95 = az.hdi(tmp, hdi_prob=0.95).mu.data
    hdi_post50 = az.hdi(tmp, hdi_prob=0.50).mu.data

    ## 描画
    fig, ax = plt.subplots(figsize=(10, 4))

    # 観測値の描画
    ax.scatter(df2['Day'], df2['Weight'], color='blue', s=30, label='観測値')

    ## 観測期間
    # 事後平均の描画
    ax.plot(df2['Day'], mean_post[:len(df2)], color='tab:blue', label='事後平均')
    # 95%HDI区間の描画
    ax.fill_between(df2['Day'], hdi_post95[:len(df2), 0], hdi_post95[:len(df2), 1],
                    color='tab:blue', alpha=0.15, label='HDI 95%')
    # 50%HDI区間の描画
    ax.fill_between(df2['Day'], hdi_post50[:len(df2), 0], hdi_post50[:len(df2), 1],
                    color='tab:blue', alpha=0.4, label='HDI 50%')

    ## 予測期間
    # 予測平均の描画
    ax.plot(idx_fore[len(df2)-1:]+1, mean_post[len(df2)-1:], color='red',
            label='予測平均')
    ax.scatter([7, 14], [mean_post[6], mean_post[13]], marker='*', s=200,
               color='red')
    # 95%HDI区間の描画：予測期間
    ax.fill_between(idx_fore[len(df2)-1:]+1, hdi_post95[len(df2)-1:, 0],
                    hdi_post95[len(df2)-1:, 1], color='tomato', alpha=0.15,
                    label='HDI 95%')
    # 50%HDI区間の描画：予測期間
    ax.fill_between(idx_fore[len(df2)-1:]+1, hdi_post50[len(df2)-1:, 0],
                    hdi_post50[len(df2)-1:, 1], color='tomato', alpha=0.4,
                    label='HDI 50%')

    # 修飾
    ax.set(xlabel='Day', ylabel='μの事後分布(kg)')
    ax.grid()
    ax.legend(bbox_to_anchor=(1, 1))
    plt.show()

描画処理を実行します。

### muの将来予測値の平均値、50%・95%HDI区間の描画 ※図4.9に相当
plot_mu_posterior2(idata1_2)

【実行結果】
テキストと近似しているように感じます（ホッとしました）。
ローカルレベルモデルはトレンド成分（傾き）や季節変動成分（周期性）を含まないので、予測の平均値は横ばい（切片的）になっている感じです。
なお赤い星印は、７日・１４日の欠損値を補った事後平均値です。
MCMC（事後分布からのサンプリング）で欠損値補間や将来予測をできることが分かりました！

⑥ 事後予測
目的変数の事後予測を行います。
テキストには無い処理です。

### 事後予測のサンプリングの実施
with model1_2:
    idata1_2.extend(pm.sample_posterior_predictive(idata1_2))

【処理結果】

目的変数「Weight」の事後予測サンプリングデータを描画します。
描画処理（予測用）を関数化します。

### 目的変数の予測値の平均値、50%・95%HDI区間の描画関数の定義

def plot_post_predict2(idata):
    
    ## 事後平均・95%HDI・50%HDIの算出
    tmp = idata.posterior_predictive
    mean_post = tmp.likelihood.stack(sample=('chain', 'draw')).mean(axis=1).data
    hdi_post95 = az.hdi(tmp, hdi_prob=0.95).likelihood.data
    hdi_post50 = az.hdi(tmp, hdi_prob=0.50).likelihood.data

    ## 描画
    fig, ax = plt.subplots(figsize=(10, 4))

    # 観測値の描画
    ax.scatter(df2['Day'], df2['Weight'], color='blue', s=30, label='観測値')

    ## 観測期間
    # 事後平均の描画
    ax.plot(df2['Day'], mean_post[:len(df2)], color='tab:blue', label='事後平均')
    # 95%HDI区間の描画
    ax.fill_between(df2['Day'], hdi_post95[:len(df2), 0], hdi_post95[:len(df2), 1],
                    color='tab:blue', alpha=0.15, label='HDI 95%')
    # 50%HDI区間の描画
    ax.fill_between(df2['Day'], hdi_post50[:len(df2), 0], hdi_post50[:len(df2), 1],
                    color='tab:blue', alpha=0.4, label='HDI 50%')

    ## 予測期間
    # 予測平均の描画
    ax.plot(idx_fore[len(df2)-1:]+1, mean_post[len(df2)-1:], color='red',
            label='予測平均')
    ax.scatter([7, 14], [mean_post[6], mean_post[13]], marker='*', s=200,
               color='red')
    # 95%HDI区間の描画：予測期間
    ax.fill_between(idx_fore[len(df2)-1:]+1, hdi_post95[len(df2)-1:, 0],
                    hdi_post95[len(df2)-1:, 1], color='tomato', alpha=0.15,
                    label='HDI 95%')
    # 50%HDI区間の描画：予測期間
    ax.fill_between(idx_fore[len(df2)-1:]+1, hdi_post50[len(df2)-1:, 0],
                    hdi_post50[len(df2)-1:, 1], color='tomato', alpha=0.4,
                    label='HDI 50%')

    # 修飾
    ax.set(xlabel='Day', ylabel='Yの事後予測(kg)')
    ax.grid()
    ax.legend(bbox_to_anchor=(1, 1))
    plt.show()

描画処理を実行します。

### 目的変数の予測値の平均値、50%・95%HDI区間の描画
plot_post_predict2(idata1_2)

【実行結果】

締めくくりは将来予測のKDEプロットの描画です。
テキストの図4.7に相当します。
plot_posteriorを利用します。
coordsで予測期間を指定することで、予測期間だけをプロットできます。

### 予測データのKDE 図4.7に相当
pm.plot_posterior(idata1_2, hdi_prob=0.95, group='posterior_predictive',
                  var_names=['likelihood'], coords={'data': range(20, 30)})
plt.tight_layout();

【実行結果】
各グラフの上部見出し２行目の数字に１を足すと予測期間になります。
例えば、右上のグラフ「20」は「21」期の予測です。

最初の１セットが終わりました。
お疲れ様です。少し休憩をとりましょう。

室内でリラックスする人のイラスト（家族）：「いらすとや」さんより


自己回帰AR(1)

状態方程式にpymcのARを用いるモデルです。
ランダムウォークは、次数１・切片$${\rho_0=0}$$・自己回帰係数$${\rho_1=1}$$の自己回帰モデル「AR(1)過程」です。
次のリンクでPyMC公式サイトのARのページを開けます。

１．基本モデル

① モデルの定義
状態「mu」 に pm.AR()を用いています。
「constant=False」で切片$${\rho_0}$$が無いことを指定します。
「rho=1」で次数分の自己回帰係数を指定します。
今回は次数１ですので係数は１個であり、自己回帰係数の値$${\rho_1=1}$$を指定します。

### モデルの定義

# モデルの定義
with pm.Model() as model2_1:
    
    ### データ関連定義
    # coordの定義
    model2_1.add_coord('data', values=df.index, mutable=True)
    # dataの定義
    y = pm.MutableData('y', value=df['Weight'].values, dims='data')
    
    ### 事前分布
    sigmaU = pm.Uniform('sigmaU', lower=0, upper=100)
    sigmaV = pm.Uniform('sigmaV', lower=0, upper=100)
    
    ### muの計算
    init_dist = pm.Normal.dist(mu=0, sigma=100)
    mu = pm.AR('mu', rho=1, sigma=sigmaU, constant=False, init_dist=init_dist,
               ar_order=1, dims='data')
    
    ### 尤度
    likelihood = pm.Normal('likelihood', mu=mu, sigma=sigmaV, observed=y,
                           dims='data')

【実行結果】出力なし

② モデルの表示・可視化

### モデルの表示
model2_1

【実行結果】
ARの引数には、自己回帰係数、標準偏差、初期分布が表示されています。

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

【実行結果】
シンプルなモデルです。

③ 事後分布からのサンプリング
MCMCを実行します。
NUTSサンプラーにnumpyroを使用するので爆速です。
実行時間は約10秒でした。

### 事後分布からのサンプリング ※NUTSサンプラーにnumpyroを使用　10秒
with model2_1:
    idata2_1 = pm.sample(draws=2000, tune=2000, chains=4, target_accept=0.9,
                         nuts_sampler='numpyro', random_seed=1234)

【実行結果（一部分）】

④ サンプリングデータの確認
$${\hat{R}}$$、事後分布の要約統計量、トレースプロットを確認します。
事後分布の収束確認はテキストにならって$${\hat{R} \leq 1.1}$$としています。

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

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

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

パラメータ等の事後分布の要約統計量です。
テキストの表4.1に相当します。

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

【実行結果】
状態「mu」の時間インデックス 0～19（時間 1～20 に相当）、および、状態ノイズsigmaU、観測ノイズsigmaVの要約統計量を表示しています。
テキストの表4.1の平均：mean、標準偏差sdの値に近いのではないでしょうか。
ただし、ガウスランダムウォークの推論値と完全には一致しないことに留意します。
なお、2.5%・97.5%に関して、この図は95%HDI区間を示しています。
一方でテキストは確信区間（信用区間）を示しており、計算基礎が異なっています。

トレースプロットを確認しましょう。

### トレースプロットの表示
pm.plot_trace(idata2_1, combined=True, figsize=(12, 8))
plt.tight_layout();

【実行結果】
収束している感じがします。
ガウスランダムウォークとは異なって、発散のバーコードが見られなくなりました。

⑤ サンプリングデータを描画
状態「mu」の事後平均値と95%HDIを描画します。
テキストの図4.4に相当します。

### muの事後分布の平均値、50%・95%HDI区間の描画 ※図4.4に相当
plot_mu_posterior(idata2_1)

【実行結果】
ガウスランダムウォークと比べて、若干ですが、95%HDIの帯が下振れ（下方移動）している感じがします。

続いて、plot_traceを用いて、sigmaUの事後分布と収束の様子を描画します。
テキストの図4.5と図4.6に相当します。

### sigma_Uの事後分布と収束の様子の描画 図4.5、4.6に相当
pm.plot_trace(idata2_1, var_names=['sigmaU'], combined=False, compact=False, 
              figsize=(10, 3), legend=True);

【実行結果】

⑥ 事後予測
目的変数「Weight」の事後予測を行います。
テキストには無い処理です。

### 事後予測のサンプリングの実施
with model2_1:
    idata2_1.extend(pm.sample_posterior_predictive(idata2_1))

【処理結果】省略

目的変数「Weight」の事後予測サンプリングデータを描画します。

### 目的変数の予測値の平均値、50%・95%HDI区間の描画
plot_post_predict(idata2_1)

【実行結果】
ガウスランダムウォークとほぼ同じ結果です。

２．将来予測モデル
欠損値の補間と１０期先予測を実施します。
テキストの図4.7、図4.9に相当する図を作成します。
将来予測モデル用に生成したデータ y_fore, idx_fore を用います。

① モデルの定義
欠損値を含むデータをモデル内の「data」で定義できません。
観測方程式（尤度）を定義する正規分布の引数 observed で観測データ y_fore を渡します。

### モデルの定義
# ※欠損値のデータはdataを利用できない。observedに目的変数を直接代入する

# モデルの定義
with pm.Model() as model2_2:
    
    ### データ関連定義
    # coordの定義
    model2_2.add_coord('data', values=idx_fore, mutable=True)
    # # dataの定義
    # y = pm.MutableData('y', value=df['Weight'].values, dims='data')
    
    ### 事前分布
    sigmaU = pm.Uniform('sigmaU', lower=0, upper=100)
    sigmaV = pm.Uniform('sigmaV', lower=0, upper=100)
    
    ### muの計算
    init_dist = pm.Normal.dist(mu=0, sigma=100)
    mu = pm.AR('mu', rho=1, sigma=sigmaU, constant=False, init_dist=init_dist,
               ar_order=1, dims='data')
    
    ### 尤度
    likelihood = pm.Normal('likelihood', mu=mu, sigma=sigmaV, observed=y_fore,
                           dims='data')

【実行結果】出力なし

② モデルの可視化

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

【実行結果】
12個の欠損値を推論する likelihood_unobserved が追加されました。

③ 事後分布からのサンプリング
MCMCを実行します。
NUTSサンプラーにnumpyroを指定するとエラーになりました。
欠損値を含むデータをnumpyroで利用する方法が分からないので、PyMC標準のサンプラーを用いることにします。
実行時間は約５分でした。

### 事後分布からのサンプリング ※NUTSサンプラーにnumpyroを未使用　5分
# 欠損値のある場合のnumpyroの使い方がわからない

with model2_2:
    idata2_2 = pm.sample(draws=2000, tune=2000, chains=4, target_accept=0.9,
                         random_seed=1234)

【実行結果】省略

④ サンプリングデータの確認
$${\hat{R}}$$、事後分布の要約統計量、トレースプロットを確認します。
事後分布の収束確認はテキストにならって$${\hat{R} \leq 1.1}$$としています。

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

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

【実行結果】
$${\hat{R}>1.1}$$のパラメータは「0」件です。
全てのパラメータが$${\hat{R} \leq 1.1}$$であることを確認できました。
また「likelihood_unobserved」および「likelihood」の推論値が追加されました。

パラメータ等の事後分布の要約統計量です。
テキストの表4.1に相当します。

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

【実行結果】
雰囲気をお楽しみ下さい。

トレースプロットを確認しましょう。

### トレースプロットの表示
pm.plot_trace(idata2_2, combined=True, figsize=(12, 8))
plt.tight_layout();

【実行結果】
収束している感じがします。
欠損値を含むこちらのサンプリングデータには発散のバーコードが見られます。

⑤ サンプリングデータを描画
将来期間を含む状態「mu」の事後平均値と95%HDIを描画します。
テキストの図4.9に相当します。

### muの将来予測値の平均値、50%・95%HDI区間の描画 ※図4.9に相当
plot_mu_posterior2(idata2_2)

【実行結果】
テキストと近似しているように感じます（ホッとしました）。
ローカルレベルモデルはトレンド成分（傾き）や季節変動成分（周期性）を含まないので、予測の平均値は横ばい（切片的）になっている感じです。
なお赤い星印は、７日・１４日の欠損値を補った事後平均値です。
MCMC（事後分布からのサンプリング）で欠損値補間や将来予測をできることが分かりました！

⑥ 事後予測
目的変数の事後予測を行います。
テキストには無い処理です。

### 事後予測のサンプリングの実施
with model2_2:
    idata2_2.extend(pm.sample_posterior_predictive(idata2_2))

【処理結果】省略

目的変数「Weight」の事後予測サンプリングデータを描画します。

### 目的変数の予測値の平均値、50%・95%HDI区間の描画
plot_post_predict2(idata2_2)

【実行結果】

締めくくりは将来予測のKDEプロットの描画です。
テキストの図4.7に相当します。
plot_posteriorを利用します。
coordsで予測期間を指定することで、予測期間だけをプロットできます。

### 予測データのKDE 図4.7に相当
pm.plot_posterior(idata2_2, hdi_prob=0.95, group='posterior_predictive',
                  var_names=['likelihood'], coords={'data': range(20, 30)})
plt.tight_layout();

【実行結果】
各グラフの上部見出し２行目の数字に１を足すと予測期間になります。
例えば、右上のグラフ「20」は「21」期の予測です。

２セット目が終わりました。
お疲れ様です。少し休憩をとりましょう。

ハイキングでお弁当を食べる家族のイラスト：「いらすとや」さんより


状態空間モデル

pymc_experimental.statespaceを用いるモデルです。
pymc_experimentalは、正式なpymc機能になる前の「実験段階のモジュール」を扱うパッケージです。
次のリンクでpymc_experimental公式サイトの状態空間モデルのページを開けます。


１．基本モデル

① モデルの定義
２段階でモデルを定義します。
１段階目は状態空間モデルの定義です。
pymc_experimental.statespaceを用いています。

### 状態空間モデルの定義

# レベル・トレンド成分の定義
# μt = μt-1 + ut,  ut ~ Normal(0, σu²)
trend = st.LevelTrendComponent(order=1, innovations_order=1)
# 測定ノイズの定義　分散パラメータσv²
error = st.MeasurementError(name='obs')
# レベル・トレンド成分と観測ノイズで状態空間モデルを構築
ss_mod = (trend + error).build(name='weight')

【補足】
■ LevelTrendComponent
ローカルレベルモデルは、状態方程式が「レベル成分」に相当します。
LevelTrendComponent＝レベル・トレンド成分で状態方程式を定義します。
引数のorder、innovations_orderの意味合いはややこしいのですが、公式サイトのガウスランダムウォーク事例に従って、それぞれ１を設定します。
時系列変数が１個かつ１次であること（order=1）、状態のノイズが１個であること（innovations_order=1）のようです。

■ MeasurementError
観測方程式の観測ノイズです。

■ (trend + error).build()
前半のカッコ内は、状態方程式を定義した trend と観測ノイズを定義した error を足して、観測方程式を完成させています。
build()は、状態空間モデルを構築するおまじないです。

【実行結果】
PyMCモデルの定義で設定の必要な「パラメータ」を指示してくれます。
指示通りのパラメータ名を用いて、形状、制約（正の数など）、次元を設定します。

では２段階目のPyMCモデルを定義しましょう。

### PyMCモデルの定義
with pm.Model(coords=ss_mod.coords) as model3_1:
    
    ## パラメータの事前分布
    # 初期値の標準偏差
    P0 = pm.Uniform('P0', lower=0, upper=1, dims=['state', 'state_aux'])
    # レベル・トレンド成分の初期値
    initial_trend = pm.Normal('initial_trend', mu=0, sigma=100,
                               dims=['trend_state'])
    # 状態変数の標準偏差
    sigma_trend = pm.Uniform('sigma_trend', lower=0, upper=100,
                             dims=['trend_shock'])
    sigma_obs = pm.Uniform('sigma_obs', lower=0, upper=100)

    ## 状態空間モデルの計算グラフ構築
    ss_mod.build_statespace_graph(data=df['Weight'], mode='JAX')

【実行結果】出力なし

【補足】
P0、initial_trend、sigma_trend、sigma_obsは状態空間モデルの定義時に指示のあったパラメータです。

  • P0：（おそらく）何かの初期値の標準偏差に使われる模様です（推測）。
    適切な設定値（お作法）がわからないので、ひとまず一様分布にしました。

  • initial_trend：レベル・トレンド成分の初期値と思われます。
    ガウスランダムウォークや自己回帰モデルの状態「mu」の初期分布と同じ内容にしました。

  • sigma_trend：状態のノイズ（システムノイズ）です。ガウスランダムウォークや自己回帰モデルのsigmaUの事前分布と同じ内容にしました。

  • sigma_obs：観測ノイズです。ガウスランダムウォークや自己回帰モデルのsigmaVの事前分布と同じ内容にしました。

「状態空間モデルの計算グラフ構築」の「ss_mod.build_statespace_graph」では、PyMCモデルの諸パラメータを使って状態空間モデル ss_mod の状態空間モデルの計算グラフ（pymcの内部構造）を生成する、こんな感じのことをしていると想像いたします。

② モデルの表示・可視化
パラメータ（変数）の名前に「_」（アンダースコア）を含むため、自環境ではTeXのエラーが発生してしまいます。
そこで、__getstate__()で内部構造を赤裸々にします！

### モデルの表示（）
model3_1.__getstate__()

【実行結果】
途中省略されたモデル内部の情報です。
読み解けません（泣）

図による可視化ならば読み解きしやすいかもです。

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

【実行結果】
読み解けません（泣）
四角ボックスの決定論的変数がおそらく、数学的な状態空間モデルに現れる各種変数と思われます。

ちなみに次の部分が「状態」です。
たぶん、状態空間モデルの「予測」predicted、「フィルタ」filtered、「平滑化」smoothedだと思われます。

③ 事後分布からのサンプリング
ここからは通常のPyMCモデルと同様の処理になります。
ではMCMCを実行しましょう！
NUTSサンプラーにnumpyroを使用するので爆速です（のはずです）。
実行時間は約55秒でした。

### 事後分布からのサンプリング ※NUTSサンプラーにnumpyroを使用　55秒
with model3_1:
    idata3_1 = pm.sample(draws=2000, tune=2000, chains=4, target_accept=0.9,
                         nuts_sampler='numpyro', random_seed=1234)

【実行結果（一部分）】

④ サンプリングデータの確認
$${\hat{R}}$$、事後分布の要約統計量、トレースプロットを確認します。
事後分布の収束確認はテキストにならって$${\hat{R} \leq 1.1}$$としています。

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

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

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

推論データ idata の外観を見てみましょう。

# 推論データの外観
idata3_1.posterior

【実行結果】
モデル等で明示的に指定していない、様々な次元（Dimensions）、様々な座標（Coordinates）、変数（Data variables）が存在します。

平滑化後の「状態」smoothed_state、状態のノイズ sigma_trend、観測ノイズ sigma_obs の事後分布の要約統計量です。
テキストの表4.1に相当します。

### 推論データの要約統計情報の表示：※表4.1に相当
var_names = ['smoothed_state', 'sigma_trend', 'sigma_obs']
pm.summary(idata3_1, hdi_prob=0.95, var_names=var_names)

【実行結果】
テキストの表4.1の平均：mean、標準偏差sdと比べてみて、状態のノイズ sigma_trend=sigma_uと観測ノイズ sigma_obs=sigma_v は相違しています。
なお、2.5%・97.5%に関して、この図は95%HDI区間を示しています。
一方でテキストは確信区間（信用区間）を示しており、計算基礎が異なっています。

トレースプロットを確認しましょう。
PyMCモデル定義で追加した４つのパラメータを見ます。

### トレースプロットの表示
pm.plot_trace(idata3_1, combined=True, var_names=ss_mod.param_names,
              figsize=(12, 10))
plt.tight_layout();

【実行結果】
収束している感じがします。
P0は一様分布の０～１の区間を幅広く漂っており、何となく良い結果とは言えない予感がします。

⑤ サンプリングデータを描画
状態smoothed_stateの事後平均値と95%HDIを描画してみます。
正直なところ、このsmoothed_stateが状態「mu」を示しているのか不安ですが、テキストの図4.4に相当するのかも、って思っています。

### 状態smoothed_stateの事後分布の平均値、50%・95%HDI区間の描画

# 事後平均・95%HDIの算出
tmp = idata3_1.posterior
mean_post = (tmp.smoothed_state.stack(sample=('chain', 'draw'))
             .squeeze().mean(axis=1))
hdi_post95 = az.hdi(tmp, hdi_prob=0.95).smoothed_state.squeeze()
hdi_post50 = az.hdi(tmp, hdi_prob=0.50).smoothed_state.squeeze()

# 描画
fig, ax = plt.subplots(figsize=(10, 4))
# 観測値の描画
ax.scatter(df['Day'], df['Weight'], color='blue', s=30, label='観測値')
# 事後平均の描画
ax.plot(df['Day'], mean_post, color='tab:blue', label='事後平均')
# 95%HDI区間の描画
ax.fill_between(df['Day'], hdi_post95[:, 0], hdi_post95[:, 1], color='tab:blue',
                alpha=0.15, label='HDI 95%')
# 50%HDI区間の描画
ax.fill_between(df['Day'], hdi_post50[:, 0], hdi_post50[:, 1], color='tab:blue',
                alpha=0.4, label='HDI 50%')
# 修飾
ax.set(xlabel='Day', ylabel='Y (kg)')
ax.grid()
ax.legend()
plt.show()

【実行結果】
なかなか尖った結果です。

続いて、plot_traceを用いて、sigma_trendの事後分布と収束の様子を描画します。
テキストの図4.5と図4.6に相当します。

### trendコンポーネントの標準偏差のトレースプロット ※図4.6に相当
pm.plot_trace(idata3_1, compact=False, var_names='sigma_trend', figsize=(10, 3))
plt.tight_layout();

【実行結果】

⑥ 状態変数の事後分布
sample_conditional_posterior()を用いて、状態変数の事後分布のサンプリングを実行します。

### 状態変数の事後分布のサンプリングの実行
post3_1 = ss_mod.sample_conditional_posterior(idata3_1, random_seed=1234)

【実行結果】省略

目的変数「Weight」の事後分布サンプリングデータを描画します。

### 目的変数の予測値の平均値、50%・95%HDI区間の描画

# 事後平均・95%HDIの算出
component_idata = ss_mod.extract_components_from_idata(post3_1)
tmp = component_idata.smoothed_posterior_observed
mean_post = tmp.stack(sample=('chain', 'draw')).squeeze().mean(axis=1)
hdi_post95 = az.hdi(tmp, hdi_prob=0.95).smoothed_posterior_observed.squeeze()
hdi_post50 = az.hdi(tmp, hdi_prob=0.50).smoothed_posterior_observed.squeeze()

# 描画
fig, ax = plt.subplots(figsize=(10, 4))
# 観測値の描画
ax.scatter(df['Day'], df['Weight'], color='blue', s=30, label='観測値')
# 事後平均の描画
ax.plot(df['Day'], mean_post, color='tab:blue', label='事後平均')
# 95%HDI区間の描画
ax.fill_between(df['Day'], hdi_post95[:, 0], hdi_post95[:, 1], color='tab:blue',
                alpha=0.15, label='HDI 95%')
# 50%HDI区間の描画
ax.fill_between(df['Day'], hdi_post50[:, 0], hdi_post50[:, 1], color='tab:blue',
                alpha=0.4, label='HDI 50%')
# 修飾
ax.set(xlabel='Day', ylabel='Y (kg)')
ax.grid()
ax.legend()
plt.show()

【実行結果】
95%HDIの帯の幅が広い感じがします。

⑦ 将来予測
pymc_experimental.statespaceのforecastを用いて将来期間の予測を実行します。
将来期間の予測が簡単に行えるのはGoodです！

### 将来予測のサンプリングの実行

# 初期値設定：予測期間
n_periods = 11

# 将来予測のサンプリングの実行
forecasts3_1 = ss_mod.forecast(idata3_1, start=df.index[-1], periods=n_periods,
                               random_seed=1234)

【実行結果】

将来予測データの外観を見てみましょう。

### 将来予測データの外観
forecasts3_1

【実行結果】
forecast_observedが目的変数の将来予測値です。

目的変数「Weight」の予測データを描画します。

### 目的変数の予測値の平均値、50%・95%HDI区間の描画

### 事後予測・将来予測の平均・HDIの算出
# 事後予測
component_idata3_1 = ss_mod.extract_components_from_idata(post3_1)
tmp = component_idata3_1.smoothed_posterior_observed
mean_post = tmp.stack(sample=('chain', 'draw')).squeeze().mean(axis=1)
hdi_post95 = az.hdi(tmp, hdi_prob=0.95).smoothed_posterior_observed.squeeze()
hdi_post50 = az.hdi(tmp, hdi_prob=0.50).smoothed_posterior_observed.squeeze()
# 将来予測
component_idata2_f = ss_mod.extract_components_from_idata(forecasts3_1)
tmp_f = component_idata2_f.forecast_observed
mean_fore = tmp_f.stack(sample=('chain', 'draw')).squeeze().mean(axis=1)
hdi_fore95 = az.hdi(tmp_f, hdi_prob=0.95).forecast_observed.squeeze()
hdi_fore50 = az.hdi(tmp_f, hdi_prob=0.50).forecast_observed.squeeze()

### 描画
fig, ax = plt.subplots(figsize=(10, 4))

## 観測値の描画
ax.scatter(df['Day'], df['Weight'], color='blue', s=30, label='観測値')

## 事後予測の描画
# 事後平均の描画
ax.plot(df['Day'], mean_post, color='tab:blue', label='事後平均')
# 95%HDI区間の描画
ax.fill_between(df['Day'], hdi_post95[:, 0], hdi_post95[:, 1], color='tab:blue',
                alpha=0.15, label='HDI 95%')
# 50%HDI区間の描画
ax.fill_between(df['Day'], hdi_post50[:, 0], hdi_post50[:, 1], color='tab:blue',
                alpha=0.4, label='HDI 50%')

## 将来予測の描画
# x軸の作成
x_add = list(range(df['Day'].max(), df['Day'].max()+11))
# 平均の描画
ax.plot(x_add, mean_fore, color='red', label='将来予測')
# 95%HDI区間の描画
ax.fill_between(x_add, hdi_fore95[:, 0], hdi_fore95[:, 1], color='tomato',
                alpha=0.15, label='HDI 95%')
# 50%HDI区間の描画
ax.fill_between(x_add, hdi_fore50[:, 0], hdi_fore50[:, 1], color='tomato',
                alpha=0.4, label='HDI 50%')

## 修飾
ax.set(xlabel='Day', ylabel='Y(kg)')
ax.grid()
ax.legend(bbox_to_anchor=(1.17, 1))
plt.show()

【実行結果】

締めくくりは将来予測のKDEプロットの描画です。
テキストの図4.7に相当します。
plot_posteriorを利用します。

### 将来予測データのKDE　20が1期先、29が20期先 ※図4.7に相当
pm.plot_posterior(forecasts3_1, hdi_prob=0.95, var_names=['forecast_observed'])
plt.tight_layout();

【実行結果】
各グラフの上部見出し２行目の数字に１を足すと予測期間になります。
例えば、上真ん中のグラフ「20」は「21」期の予測です。
実績期間の最終日である20日（Pythonのインデックス上は19）も将来予測に含めているので、19のグラフが最初に描画されています。

２．将来予測
欠損値の補間と１０期先予測を実施します。
テキストの図4.7、図4.9に相当する図を作成します。
将来予測モデル用に生成したデータ df2 を用います。
将来期間の予測機能があるので、将来期間を欠損値で埋めたデータを設定しなくてもＯＫです。

① モデルの定義
２段階でモデルを定義します。
１段階目は状態空間モデルの定義です。
pymc_experimental.statespaceを用います。
ただし、「１．基本モデル」の定義をそのまま使えるので、改めてモデルをしないことにします。

では２段階目のPyMCモデルを定義しましょう。
「１．基本モデル」のモデルと区別しておきたいので、モデル定義を行います。
「１．基本モデル」と異なるのは最後の１文、観測データを指定する「data=df2['Weight']」です。

### PyMCモデルの定義
with pm.Model(coords=ss_mod.coords) as model3_2:
    
    ## パラメータの事前分布
    # 初期値の標準偏差
    P0 = pm.Uniform('P0', lower=0, upper=1, dims=['state', 'state_aux'])
    # レベル・トレンド成分の初期値
    initial_trend = pm.Normal('initial_trend', mu=0, sigma=100,
                              dims=['trend_state'])
    # 状態変数の標準偏差
    sigma_trend = pm.Uniform('sigma_trend', lower=0, upper=100,
                             dims=['trend_shock'])
    sigma_obs = pm.Uniform('sigma_obs', lower=0, upper=100)

    ## 状態空間モデルの計算グラフ構築
    ss_mod.build_statespace_graph(data=df2['Weight'], mode='JAX')

【実行結果】出力なし

② モデルの表示・可視化
「１．基本モデル」と大差ないと思われますので、実行を省略します。

③ 事後分布からのサンプリング
ではMCMCを実行しましょう！
NUTSサンプラーにnumpyroを使用するので爆速です（のはずです）。
実行時間は約50秒でした。

### 事後分布からのサンプリング ※NUTSサンプラーにnumpyroを使用　50秒
# 欠損値が含まれていてもnumpyroを利用できる

with model3_2:
    idata3_2 = pm.sample(draws=2000, tune=2000, chains=4, target_accept=0.9,
                         nuts_sampler='numpyro', random_seed=1234)

【実行結果（一部分）】

④ サンプリングデータの確認
$${\hat{R}}$$、事後分布の要約統計量、トレースプロットを確認します。
事後分布の収束確認はテキストにならって$${\hat{R} \leq 1.1}$$としています。

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

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

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

平滑化後の「状態」smoothed_state、状態のノイズ sigma_trend、観測ノイズ sigma_obs の事後分布の要約統計量です。
テキストの表4.1に相当します。

### 推論データの要約統計情報の表示：※表4.1に相当
var_names = ['smoothed_state', 'sigma_trend', 'sigma_obs']
pm.summary(idata3_2, hdi_prob=0.95, var_names=var_names)

【実行結果】
観測データが欠損値となっている７日（インデックス６）、１４日（インデックス１３）にも「状態」の推論値が設定されています。

トレースプロットを確認しましょう。
PyMCモデル定義で追加した４つのパラメータを見ます。

### トレースプロットの表示
pm.plot_trace(idata3_2, combined=True, var_names=ss_mod.param_names,
              figsize=(12, 10))
plt.tight_layout();

【実行結果】
収束している感じがします。
P0は一様分布の０～１の区間を幅広く漂っており、何となく良い結果とは言えない予感がします。
観測データに欠損値を含んでいることが原因かもですが、発散を示すバーコードが見られます。

⑤ サンプリングデータを描画
状態smoothed_stateの事後平均値と95%HDIを描画してみます。
正直なところ、このsmoothed_stateが状態「mu」を示しているのか不安ですが、テキストの図4.4に相当するのかも、って思っています。

### 状態smoothed_stateの事後分布の平均値、50%・95%HDI区間の描画

# 事後平均・95%HDIの算出
tmp = idata3_2.posterior
mean_post = (tmp.smoothed_state.stack(sample=('chain', 'draw'))
             .squeeze().mean(axis=1))
hdi_post95 = az.hdi(tmp, hdi_prob=0.95).smoothed_state.squeeze()
hdi_post50 = az.hdi(tmp, hdi_prob=0.50).smoothed_state.squeeze()

# 描画
fig, ax = plt.subplots(figsize=(10, 4))
# 観測値の描画
ax.scatter(df2['Day'], df2['Weight'], color='blue', s=30, label='観測値')
# 事後平均の描画
ax.plot(df2['Day'], mean_post, color='tab:blue', label='事後平均')
ax.scatter([7, 14], [mean_post[6], mean_post[13]], marker='*', s=200, color='red')
# 95%HDI区間の描画
ax.fill_between(df2['Day'], hdi_post95[:, 0], hdi_post95[:, 1], color='tab:blue',
                alpha=0.15, label='HDI 95%')
# 50%HDI区間の描画
ax.fill_between(df2['Day'], hdi_post50[:, 0], hdi_post50[:, 1], color='tab:blue',
                alpha=0.4, label='HDI 50%')
# 修飾
ax.set(xlabel='Day', ylabel='Y (kg)')
ax.grid()
ax.legend()
plt.show()

【実行結果】
「１．基本モデル」と比べて、欠損値部分の95%HDIの幅が広い感じがします（特に14日目）。

続いて、plot_traceを用いて、sigma_trendの事後分布と収束の様子を描画します。
テキストの図4.5と図4.6に相当します。

### trendコンポーネントの標準偏差のトレースプロット ※図4.5、4.6に相当
pm.plot_trace(idata3_2, compact=False, var_names='sigma_trend', figsize=(10, 3))
plt.tight_layout();

【実行結果】

⑥ 状態変数の事後分布
sample_conditional_posterior()を用いて、状態変数の事後分布のサンプリングを実行します。

### 状態変数の事後分布のサンプリングの実行
post3_2 = ss_mod.sample_conditional_posterior(idata3_2, random_seed=1234)

【実行結果】省略

目的変数「Weight」の事後分布サンプリングデータを描画します。

### 目的変数の予測値の平均値、50%・95%HDI区間の描画

# 事後平均・95%HDIの算出
component_idata = ss_mod.extract_components_from_idata(post3_2)
tmp = component_idata.smoothed_posterior_observed
mean_post = tmp.stack(sample=('chain', 'draw')).squeeze().mean(axis=1)
hdi_post95 = az.hdi(tmp, hdi_prob=0.95).smoothed_posterior_observed.squeeze()
hdi_post50 = az.hdi(tmp, hdi_prob=0.50).smoothed_posterior_observed.squeeze()

# 描画
fig, ax = plt.subplots(figsize=(10, 4))
# 観測値の描画
ax.scatter(df2['Day'], df2['Weight'], color='blue', s=30, label='観測値')
# 事後平均の描画
ax.plot(df2['Day'], mean_post, color='tab:blue', label='事後平均')
ax.scatter([7, 14], [mean_post[6], mean_post[13]], marker='*', s=200, color='red')
# 95%HDI区間の描画
ax.fill_between(df2['Day'], hdi_post95[:, 0], hdi_post95[:, 1], color='tab:blue',
                alpha=0.15, label='HDI 95%')
# 50%HDI区間の描画
ax.fill_between(df2['Day'], hdi_post50[:, 0], hdi_post50[:, 1], color='tab:blue',
                alpha=0.4, label='HDI 50%')
# 修飾
ax.set(xlabel='Day', ylabel='Y (kg)')
ax.grid()
ax.legend()
plt.show()

【実行結果】
欠損値部分（赤い星）の95%HDIの帯の幅が広い感じがします。

⑦ 将来予測
pymc_experimental.statespaceのforecastを用いて将来期間の予測を実行します。
将来期間の予測が簡単に行えるのはGoodです！

### 将来予測のサンプリングの実行

# 初期値設定：予測期間
n_periods = 11

# 将来予測のサンプリングの実行
forecasts3_2 = ss_mod.forecast(idata3_2, start=df.index[-1], periods=n_periods,
                               random_seed=1234)

【実行結果】

目的変数「Weight」の予測データを描画します。

### 目的変数の予測値の平均値、50%・95%HDI区間の描画

### 事後予測・将来予測の平均・HDIの算出
# 事後予測
component_idata3_2 = ss_mod.extract_components_from_idata(post3_2)
tmp = component_idata3_2.smoothed_posterior_observed
mean_post = tmp.stack(sample=('chain', 'draw')).squeeze().mean(axis=1)
hdi_post95 = az.hdi(tmp, hdi_prob=0.95).smoothed_posterior_observed.squeeze()
hdi_post50 = az.hdi(tmp, hdi_prob=0.50).smoothed_posterior_observed.squeeze()
# 将来予測
component_idata2_f = ss_mod.extract_components_from_idata(forecasts3_2)
tmp_f = component_idata2_f.forecast_observed
mean_fore = tmp_f.stack(sample=('chain', 'draw')).squeeze().mean(axis=1)
hdi_fore95 = az.hdi(tmp_f, hdi_prob=0.95).forecast_observed.squeeze()
hdi_fore50 = az.hdi(tmp_f, hdi_prob=0.50).forecast_observed.squeeze()

### 描画
fig, ax = plt.subplots(figsize=(10, 4))

## 観測値の描画
ax.scatter(df2['Day'], df2['Weight'], color='blue', s=30, label='観測値')

## 事後予測の描画
# 事後平均の描画
ax.plot(df2['Day'], mean_post, color='tab:blue', label='事後平均')
# 95%HDI区間の描画
ax.fill_between(df2['Day'], hdi_post95[:, 0], hdi_post95[:, 1], color='tab:blue',
                alpha=0.15, label='HDI 95%')
# 50%HDI区間の描画
ax.fill_between(df2['Day'], hdi_post50[:, 0], hdi_post50[:, 1], color='tab:blue',
                alpha=0.4, label='HDI 50%')

## 将来予測の描画
# x軸の作成
x_add = list(range(df2['Day'].max(), df2['Day'].max()+11))
# 平均の描画
ax.plot(x_add, mean_fore, color='red', label='将来予測')
ax.scatter([7, 14], [mean_post[6], mean_post[13]], marker='*', s=200, color='red')
# 95%HDI区間の描画
ax.fill_between(x_add, hdi_fore95[:, 0], hdi_fore95[:, 1], color='tomato',
                alpha=0.15, label='HDI 95%')
# 50%HDI区間の描画
ax.fill_between(x_add, hdi_fore50[:, 0], hdi_fore50[:, 1], color='tomato',
                alpha=0.4, label='HDI 50%')

## 修飾
ax.set(xlabel='Day', ylabel='Y(kg)')
ax.grid()
ax.legend(bbox_to_anchor=(1.17, 1))
plt.show()

【実行結果】

締めくくりは将来予測のKDEプロットの描画です。
テキストの図4.7に相当します。
plot_posteriorを利用します。

### 将来予測データのKDE　20が1期先、29が20期先 ※図4.7に相当
pm.plot_posterior(forecasts3_2, hdi_prob=0.95, var_names=['forecast_observed'])
plt.tight_layout();

【実行結果】
各グラフの上部見出し２行目の数字に１を足すと予測期間になります。
例えば、上真ん中のグラフ「20」は「21」期の予測です。
実績期間の最終日である20日（Pythonのインデックス上は19）も将来予測に含めているので、19のグラフが最初に描画されています。

最終セットが終わりました。
お疲れ様でした。

からくり人形のイラスト：「いらすとや」さんより

まとめ

  • 状態空間モデルのシンプルなケース「ローカルレベルモデル」を検討しました。

  • ローカルレベルモデルは状態方程式が「ランダムウォーク」なので、通常のPyMCの「ガウスランダムウォーク」または「自己回帰AR」でモデリングできることが分かりました。

  • 通常のPyMCで将来期間を予測する方法の１つに「将来予測期間に見合った欠損値データを渡す」対応方法があることが分かりました。

  • 「pymc_experimental.statespace」で「ガウス状態空間時系列モデル」を構築でき、将来予測機能があることが分かりました。

そうそう、推論データ（idata）を保存いたしましょう。
pickle形式で保存します。

### idataの保存 pickle
import pickle

file = r'ch04_01_idata1_1.pkl'
with open(file, 'wb') as f:
    pickle.dump(idata1_1, f)
    
file = r'ch04_01_idata1_2.pkl'
with open(file, 'wb') as f:
    pickle.dump(idata1_2, f)

file = r'ch04_01_idata2_1.pkl'
with open(file, 'wb') as f:
    pickle.dump(idata2_1, f)
    
file = r'ch04_01_idata2_2.pkl'
with open(file, 'wb') as f:
    pickle.dump(idata2_2, f)

file = r'ch04_01_idata3_1.pkl'
with open(file, 'wb') as f:
    pickle.dump(idata3_1, f)

file = r'ch04_01_idata3_2.pkl'
with open(file, 'wb') as f:
    pickle.dump(idata3_2, f)

読み込みコードは次のとおりです。
「file」のパス＆ファイル名と「idata1_1_load」を適宜変更してお使い下さい。

### idataの読み込み（例）
file = r'ch04_01_idata1_1.pkl'
with open(file, 'rb') as f:
    idata1_1_load = pickle.load(f)


第４章 その２ へ続く

■次回の取り組みテーマ
３つの状態空間モデル
・レベルの変化のみを考慮したローカルレベルモデル
・レベルの変化に加えて傾きの変化も考慮したローカル線形トレンドモデル
・季節成分を考慮したローカルレベルモデル

次の記事

前の記事

目次


ブログの紹介

note で４つのシリーズ記事を書いています。
ぜひ覗いていってくださいね！

１．のんびり統計
統計検定２級の問題集を手がかりにして、確率・統計をざっくり掘り下げるブログです。
雑談感覚で大丈夫です。ぜひ覗いていってくださいね。
統計検定２級公式問題集CBT対応版に対応しています。

２．実験！たのしいベイズモデリングをPyMC Ver.5で

書籍「たのしいベイズモデリング」の心理学研究に用いられたベイズモデルを PyMC Ver.5で描いて分析します。
この書籍をはじめ、多くのベイズモデルは Ｒ言語＋Stanで書かれています。
PyMCの可能性を探り出し、手軽にベイズモデリングを実践できるように努めます。
身近なテーマ、イメージしやすいテーマですので、ぜひぜひPyMCで動かして、一緒に楽しみましょう！

３．Python機械学習プログラミング実践記
書籍「Python機械学習プログラミング PyTorch & scikit-learn編」を学んだときのさまざまな思いを記事にしました。
この書籍は、scikit-learn と PyTorch の教科書です。
よかったらぜひ、お試しくださいませ。

４．データサイエンスっぽいことを綴る
統計、データ分析、AI、機械学習、Python のコラムを不定期に綴っています。
「統計」「Python」「数学とPython」「R」のシリーズが生まれています。

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

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

お金について考える

35,930件

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