カスタム状態空間モデル

状態空間モデルの真の力は、カスタムモデルの作成と推定を可能にすることです。このノートブックでは、sm.tsa.statespace.MLEModel をサブクラス化するさまざまな状態空間モデルを示します。

一般的な状態空間モデルは、次の一般的な方法で記述できることを思い出してください。

\[\begin{split}\begin{aligned} y_t & = Z_t \alpha_{t} + d_t + \varepsilon_t \\ \alpha_{t+1} & = T_t \alpha_{t} + c_t + R_t \eta_{t} \end{aligned}\end{split}\]

オブジェクトの詳細と次元は、このリンク で確認できます。

ほとんどのモデルには、これらの要素がすべて含まれているわけではありません。たとえば、計画行列 \(Z_t\) は時間依存しない場合があり (\(\forall t \;Z_t = Z\))、モデルには観測切片 \(d_t\) がない場合があります。

まず比較的単純なものから始め、それを少しずつ拡張してより多くの要素を含める方法を示します。

  • モデル1：時間変動係数。2つの状態方程式を持つ1つの観測方程式

  • モデル2：非単位行列の遷移行列を持つ時間変動パラメーター

  • モデル3：複数の観測方程式と複数の状態方程式

  • ボーナス：ベイズ推定のためのpymc3

[ ]:

%matplotlib inline

from collections import OrderedDict

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import statsmodels.api as sm

plt.rc("figure", figsize=(16, 8))
plt.rc("font", size=15)

モデル1：時間変動係数

\[\begin{split}\begin{aligned} y_t & = d + x_t \beta_{x,t} + w_t \beta_{w,t} + \varepsilon_t \hspace{4em} \varepsilon_t \sim N(0, \sigma_\varepsilon^2)\\ \begin{bmatrix} \beta_{x,t} \\ \beta_{w,t} \end{bmatrix} & = \begin{bmatrix} \beta_{x,t-1} \\ \beta_{w,t-1} \end{bmatrix} + \begin{bmatrix} \zeta_{x,t} \\ \zeta_{w,t} \end{bmatrix} \hspace{3.7em} \begin{bmatrix} \zeta_{x,t} \\ \zeta_{w,t} \end{bmatrix} \sim N \left ( \begin{bmatrix} 0 \\ 0 \end{bmatrix}, \begin{bmatrix} \sigma_{\beta, x}^2 & 0 \\ 0 & \sigma_{\beta, w}^2 \end{bmatrix} \right ) \end{aligned}\end{split}\]

観測されたデータは \(y_t, x_t, w_t\) です。\(x_t, w_t\) は外生変数です。計画行列は時間変動であるため、3つの次元 (k_endog x k_states x nobs) を持つことに注意してください。

状態は \(\beta_{x,t}\) と \(\beta_{w,t}\) です。状態方程式は、これらの状態がランダムウォークで変化することを示しています。したがって、この場合、遷移行列は 2 x 2 の単位行列です。

最初にデータをシミュレーションし、モデルを構築し、最後にそれを推定します。

[ ]:

def gen_data_for_model1():
    nobs = 1000

    rs = np.random.RandomState(seed=93572)

    d = 5
    var_y = 5
    var_coeff_x = 0.01
    var_coeff_w = 0.5

    x_t = rs.uniform(size=nobs)
    w_t = rs.uniform(size=nobs)
    eps = rs.normal(scale=var_y ** 0.5, size=nobs)

    beta_x = np.cumsum(rs.normal(size=nobs, scale=var_coeff_x ** 0.5))
    beta_w = np.cumsum(rs.normal(size=nobs, scale=var_coeff_w ** 0.5))

    y_t = d + beta_x * x_t + beta_w * w_t + eps
    return y_t, x_t, w_t, beta_x, beta_w


y_t, x_t, w_t, beta_x, beta_w = gen_data_for_model1()
_ = plt.plot(y_t)

[ ]:

class TVRegression(sm.tsa.statespace.MLEModel):
    def __init__(self, y_t, x_t, w_t):
        exog = np.c_[x_t, w_t]  # shaped nobs x 2

        super(TVRegression, self).__init__(
            endog=y_t, exog=exog, k_states=2, initialization="diffuse"
        )

        # Since the design matrix is time-varying, it must be
        # shaped k_endog x k_states x nobs
        # Notice that exog.T is shaped k_states x nobs, so we
        # just need to add a new first axis with shape 1
        self.ssm["design"] = exog.T[np.newaxis, :, :]  # shaped 1 x 2 x nobs
        self.ssm["selection"] = np.eye(self.k_states)
        self.ssm["transition"] = np.eye(self.k_states)

        # Which parameters need to be positive?
        self.positive_parameters = slice(1, 4)

    @property
    def param_names(self):
        return ["intercept", "var.e", "var.x.coeff", "var.w.coeff"]

    @property
    def start_params(self):
        """
        Defines the starting values for the parameters
        The linear regression gives us reasonable starting values for the constant
        d and the variance of the epsilon error
        """
        exog = sm.add_constant(self.exog)
        res = sm.OLS(self.endog, exog).fit()
        params = np.r_[res.params[0], res.scale, 0.001, 0.001]
        return params

    def transform_params(self, unconstrained):
        """
        We constraint the last three parameters
        ('var.e', 'var.x.coeff', 'var.w.coeff') to be positive,
        because they are variances
        """
        constrained = unconstrained.copy()
        constrained[self.positive_parameters] = (
            constrained[self.positive_parameters] ** 2
        )
        return constrained

    def untransform_params(self, constrained):
        """
        Need to unstransform all the parameters you transformed
        in the `transform_params` function
        """
        unconstrained = constrained.copy()
        unconstrained[self.positive_parameters] = (
            unconstrained[self.positive_parameters] ** 0.5
        )
        return unconstrained

    def update(self, params, **kwargs):
        params = super(TVRegression, self).update(params, **kwargs)

        self["obs_intercept", 0, 0] = params[0]
        self["obs_cov", 0, 0] = params[1]
        self["state_cov"] = np.diag(params[2:4])

次に、カスタムモデルクラスで推定します

[ ]:

mod = TVRegression(y_t, x_t, w_t)
res = mod.fit()

print(res.summary())

データを生成した値は次のとおりです

  • 切片 = 5

  • var.e = 5

  • var.x.coeff = 0.01

  • var.w.coeff = 0.5

ご覧のとおり、推定値は実際のパラメーターをかなりうまく回復しています。

基礎となる係数 (またはカルマンフィルターでの状態) の推定された変化を回復することもできます。

[ ]:

fig, axes = plt.subplots(2, figsize=(16, 8))

ss = pd.DataFrame(res.smoothed_state.T, columns=["x", "w"])

axes[0].plot(beta_x, label="True")
axes[0].plot(ss["x"], label="Smoothed estimate")
axes[0].set(title="Time-varying coefficient on x_t")
axes[0].legend()

axes[1].plot(beta_w, label="True")
axes[1].plot(ss["w"], label="Smoothed estimate")
axes[1].set(title="Time-varying coefficient on w_t")
axes[1].legend()

fig.tight_layout();

モデル2：非単位行列の遷移行列を持つ時間変動パラメーター

これはモデル1からの小さな拡張です。単位遷移行列を持つ代わりに、推定する必要がある2つのパラメーター (\(\rho_1, \rho_2\)) を持つ行列を使用します。

\[\begin{split}\begin{aligned} y_t & = d + x_t \beta_{x,t} + w_t \beta_{w,t} + \varepsilon_t \hspace{4em} \varepsilon_t \sim N(0, \sigma_\varepsilon^2)\\ \begin{bmatrix} \beta_{x,t} \\ \beta_{w,t} \end{bmatrix} & = \begin{bmatrix} \rho_1 & 0 \\ 0 & \rho_2 \end{bmatrix} \begin{bmatrix} \beta_{x,t-1} \\ \beta_{w,t-1} \end{bmatrix} + \begin{bmatrix} \zeta_{x,t} \\ \zeta_{w,t} \end{bmatrix} \hspace{3.7em} \begin{bmatrix} \zeta_{x,t} \\ \zeta_{w,t} \end{bmatrix} \sim N \left ( \begin{bmatrix} 0 \\ 0 \end{bmatrix}, \begin{bmatrix} \sigma_{\beta, x}^2 & 0 \\ 0 & \sigma_{\beta, w}^2 \end{bmatrix} \right ) \end{aligned}\end{split}\]

前回のクラスで何を修正して動作させる必要がありますか？ + 良いニュース：多くはありません！ + 悪いニュース：いくつかのことに注意する必要があります

1) 開始パラメーター関数の変更

新しいパラメーター \(\rho_1, \rho_2\) の名前を追加する必要があり、対応する開始値を設定する必要があります。

param_names 関数は次のように変わります

def param_names(self):
    return ['intercept', 'var.e', 'var.x.coeff', 'var.w.coeff']

から

def param_names(self):
    return ['intercept', 'var.e', 'var.x.coeff', 'var.w.coeff',
           'rho1', 'rho2']

また、start_params 関数を次のように変更します

def start_params(self):
    exog = sm.add_constant(self.exog)
    res = sm.OLS(self.endog, exog).fit()
    params = np.r_[res.params[0], res.scale, 0.001, 0.001]
    return params

から

def start_params(self):
    exog = sm.add_constant(self.exog)
    res = sm.OLS(self.endog, exog).fit()
    params = np.r_[res.params[0], res.scale, 0.001, 0.001, 0.8, 0.8]
    return params

  1. update 関数の変更

次のように変わります

def update(self, params, **kwargs):
    params = super(TVRegression, self).update(params, **kwargs)

    self['obs_intercept', 0, 0] = params[0]
    self['obs_cov', 0, 0] = params[1]
    self['state_cov'] = np.diag(params[2:4])

から

def update(self, params, **kwargs):
    params = super(TVRegression, self).update(params, **kwargs)

    self['obs_intercept', 0, 0] = params[0]
    self['obs_cov', 0, 0] = params[1]
    self['state_cov'] = np.diag(params[2:4])
    self['transition', 0, 0] = params[4]
    self['transition', 1, 1] = params[5]

  1. (オプション) transform_params と untransform_params の変更

これは必須ではありませんが、\(\rho_1, \rho_2\) を -1 から 1 の間に制限したい場合があります。その場合、最初に statsmodels から 2 つのユーティリティ関数をインポートします。

from statsmodels.tsa.statespace.tools import (
    constrain_stationary_univariate, unconstrain_stationary_univariate)

constrain_stationary_univariate は、値を -1 から 1 の範囲内に制約します。unconstrain_stationary_univariate は逆関数を提供します。変換および非変換パラメーター関数は次のようになります (\(\rho_1, \rho_2\) は 4 番目と 5 番目のインデックスにあることに注意してください)。

def transform_params(self, unconstrained):
    constrained = unconstrained.copy()
    constrained[self.positive_parameters] = constrained[self.positive_parameters]**2
    constrained[4] = constrain_stationary_univariate(constrained[4:5])
    constrained[5] = constrain_stationary_univariate(constrained[5:6])
    return constrained

def untransform_params(self, constrained):
    unconstrained = constrained.copy()
    unconstrained[self.positive_parameters] = unconstrained[self.positive_parameters]**0.5
    unconstrained[4] = unconstrain_stationary_univariate(constrained[4:5])
    unconstrained[5] = unconstrain_stationary_univariate(constrained[5:6])
    return unconstrained

以下に完全なクラス (先ほど説明したオプションの変更なし) を記述します。

[ ]:

class TVRegressionExtended(sm.tsa.statespace.MLEModel):
    def __init__(self, y_t, x_t, w_t):
        exog = np.c_[x_t, w_t]  # shaped nobs x 2

        super(TVRegressionExtended, self).__init__(
            endog=y_t, exog=exog, k_states=2, initialization="diffuse"
        )

        # Since the design matrix is time-varying, it must be
        # shaped k_endog x k_states x nobs
        # Notice that exog.T is shaped k_states x nobs, so we
        # just need to add a new first axis with shape 1
        self.ssm["design"] = exog.T[np.newaxis, :, :]  # shaped 1 x 2 x nobs
        self.ssm["selection"] = np.eye(self.k_states)
        self.ssm["transition"] = np.eye(self.k_states)

        # Which parameters need to be positive?
        self.positive_parameters = slice(1, 4)

    @property
    def param_names(self):
        return ["intercept", "var.e", "var.x.coeff", "var.w.coeff", "rho1", "rho2"]

    @property
    def start_params(self):
        """
        Defines the starting values for the parameters
        The linear regression gives us reasonable starting values for the constant
        d and the variance of the epsilon error
        """

        exog = sm.add_constant(self.exog)
        res = sm.OLS(self.endog, exog).fit()
        params = np.r_[res.params[0], res.scale, 0.001, 0.001, 0.7, 0.8]
        return params

    def transform_params(self, unconstrained):
        """
        We constraint the last three parameters
        ('var.e', 'var.x.coeff', 'var.w.coeff') to be positive,
        because they are variances
        """
        constrained = unconstrained.copy()
        constrained[self.positive_parameters] = (
            constrained[self.positive_parameters] ** 2
        )
        return constrained

    def untransform_params(self, constrained):
        """
        Need to unstransform all the parameters you transformed
        in the `transform_params` function
        """
        unconstrained = constrained.copy()
        unconstrained[self.positive_parameters] = (
            unconstrained[self.positive_parameters] ** 0.5
        )
        return unconstrained

    def update(self, params, **kwargs):
        params = super(TVRegressionExtended, self).update(params, **kwargs)

        self["obs_intercept", 0, 0] = params[0]
        self["obs_cov", 0, 0] = params[1]
        self["state_cov"] = np.diag(params[2:4])
        self["transition", 0, 0] = params[4]
        self["transition", 1, 1] = params[5]

推定するには、モデル1と同じデータを使用し、\(\rho_1, \rho_2\) が1に近いと予想します。

結果は非常に良好に見えます！この推定は、\(\rho_1, \rho_2\) の開始値に非常に敏感である可能性があることに注意してください。より低い値を試すと、収束に失敗することがわかります。

[ ]:

mod = TVRegressionExtended(y_t, x_t, w_t)
res = mod.fit(maxiter=2000)  # it doesn't converge with 50 iters
print(res.summary())

モデル3：複数の観測方程式と状態方程式

時間変動パラメーターは維持しますが、今回は2つの観測方程式も使用します。

観測方程式

\(\hat{i_t}, \hat{M_t}, \hat{s_t}\) は各期間に観測されます。

観測方程式のモデルには2つの方程式があります

\[\hat{i_t} = \alpha_1 * \hat{s_t} + \varepsilon_1\]
\[\hat{M_t} = \alpha_2 + \varepsilon_2\]

状態空間モデルの一般的な表記 に従うと、観測方程式の内生部分は \(y_t = (\hat{i_t}, \hat{M_t})\) であり、外生変数は1つだけ \(\hat{s_t}\) です。

状態方程式

\[\alpha_{1, t+1} = \delta_1 \alpha_{1, t} + \delta_2 \alpha_{2, t} + W_1\]
\[\alpha_{2, t+1} = \delta_3 \alpha_{2, t} + W_2\]

状態空間モデルの行列表記

\[\begin{split}\begin{aligned} \begin{bmatrix} \hat{i_t} \\ \hat{M_t} \end{bmatrix} &= \begin{bmatrix} \hat{s_t} & 0 \\ 0 & 1 \end{bmatrix} \begin{bmatrix} \alpha_{1, t} \\ \alpha_{2, t} \end{bmatrix} + \begin{bmatrix} \varepsilon_{1, t} \\ \varepsilon_{1, t} \end{bmatrix} \hspace{6.5em} \varepsilon_t \sim N \left ( \begin{bmatrix} 0 \\ 0 \end{bmatrix}, \begin{bmatrix} \sigma_{\varepsilon_1}^2 & 0 \\ 0 & \sigma_{\varepsilon_2}^2 \end{bmatrix} \right ) \\ \begin{bmatrix} \alpha_{1, t+1} \\ \alpha_{2, t+1} \end{bmatrix} & = \begin{bmatrix} \delta_1 & \delta_1 \\ 0 & \delta_3 \end{bmatrix} \begin{bmatrix} \alpha_{1, t} \\ \alpha_{2, t} \end{bmatrix} + \begin{bmatrix} W_1 \\ W_2 \end{bmatrix} \hspace{3.em} \begin{bmatrix} W_1 \\ W_2 \end{bmatrix} \sim N \left ( \begin{bmatrix} 0 \\ 0 \end{bmatrix}, \begin{bmatrix} \sigma_{W_1}^2 & 0 \\ 0 & \sigma_{W_2}^2 \end{bmatrix} \right ) \end{aligned}\end{split}\]

いくつかのデータをシミュレートし、何を修正する必要があるかについて説明し、最後にモデルを推定して、何か妥当なものを回復しているかどうかを確認します。

[ ]:

true_values = {
    "var_e1": 0.01,
    "var_e2": 0.01,
    "var_w1": 0.01,
    "var_w2": 0.01,
    "delta1": 0.8,
    "delta2": 0.5,
    "delta3": 0.7,
}


def gen_data_for_model3():
    # Starting values
    alpha1_0 = 2.1
    alpha2_0 = 1.1

    t_max = 500

    def gen_i(alpha1, s):
        return alpha1 * s + np.sqrt(true_values["var_e1"]) * np.random.randn()

    def gen_m_hat(alpha2):
        return 1 * alpha2 + np.sqrt(true_values["var_e2"]) * np.random.randn()

    def gen_alpha1(alpha1, alpha2):
        w1 = np.sqrt(true_values["var_w1"]) * np.random.randn()
        return true_values["delta1"] * alpha1 + true_values["delta2"] * alpha2 + w1

    def gen_alpha2(alpha2):
        w2 = np.sqrt(true_values["var_w2"]) * np.random.randn()
        return true_values["delta3"] * alpha2 + w2

    s_t = 0.3 + np.sqrt(1.4) * np.random.randn(t_max)
    i_hat = np.empty(t_max)
    m_hat = np.empty(t_max)

    current_alpha1 = alpha1_0
    current_alpha2 = alpha2_0
    for t in range(t_max):
        # Obs eqns
        i_hat[t] = gen_i(current_alpha1, s_t[t])
        m_hat[t] = gen_m_hat(current_alpha2)

        # state eqns
        new_alpha1 = gen_alpha1(current_alpha1, current_alpha2)
        new_alpha2 = gen_alpha2(current_alpha2)

        # Update states for next period
        current_alpha1 = new_alpha1
        current_alpha2 = new_alpha2

    return i_hat, m_hat, s_t


i_hat, m_hat, s_t = gen_data_for_model3()

何を修正する必要がありますか？

もう一度、あまり変更する必要はありませんが、次元には注意する必要があります。

1) __init__ 関数は次のように変更します

def __init__(self, y_t, x_t, w_t):
        exog = np.c_[x_t, w_t]

        super(TVRegressionExtended, self).__init__(
            endog=y_t, exog=exog, k_states=2,
            initialization='diffuse')

        self.ssm['design'] = exog.T[np.newaxis, :, :]  # shaped 1 x 2 x nobs
        self.ssm['selection'] = np.eye(self.k_states)
        self.ssm['transition'] = np.eye(self.k_states)

から

def __init__(self, i_t: np.array, s_t: np.array, m_t: np.array):

        exog = np.c_[s_t, np.repeat(1, len(s_t))]  # exog.shape => (nobs, 2)

        super(MultipleYsModel, self).__init__(
            endog=np.c_[i_t, m_t], exog=exog, k_states=2,
            initialization='diffuse')

        self.ssm['design'] = np.zeros((self.k_endog, self.k_states, self.nobs))
        self.ssm['design', 0, 0, :] = s_t
        self.ssm['design', 1, 1, :] = 1

k_endog をどこにも指定する必要がなかったことに注意してください。初期化では、endog 行列の次元を確認した後、これを行います。

2) update() 関数

次のように変更します

def update(self, params, **kwargs):
    params = super(TVRegressionExtended, self).update(params, **kwargs)

    self['obs_intercept', 0, 0] = params[0]
    self['obs_cov', 0, 0] = params[1]

    self['state_cov'] = np.diag(params[2:4])
    self['transition', 0, 0] = params[4]
    self['transition', 1, 1] = params[5]

から

def update(self, params, **kwargs):
    params = super(MultipleYsModel, self).update(params, **kwargs)


    #The following line is not needed (by default, this matrix is initialized by zeroes),
    #But I leave it here so the dimensions are clearer
    self['obs_intercept'] = np.repeat([np.array([0, 0])], self.nobs, axis=0).T
    self['obs_cov', 0, 0] = params[0]
    self['obs_cov', 1, 1] = params[1]

    self['state_cov'] = np.diag(params[2:4])
    #delta1, delta2, delta3
    self['transition', 0, 0] = params[4]
    self['transition', 0, 1] = params[5]
    self['transition', 1, 1] = params[6]

残りのメソッドは、非常に明白な方法で変更されます (パラメーター名を追加し、インデックスが機能することを確認するなど)。関数の完全なコードは次のとおりです。

[ ]:

starting_values = {
    "var_e1": 0.2,
    "var_e2": 0.1,
    "var_w1": 0.15,
    "var_w2": 0.18,
    "delta1": 0.7,
    "delta2": 0.1,
    "delta3": 0.85,
}


class MultipleYsModel(sm.tsa.statespace.MLEModel):
    def __init__(self, i_t: np.array, s_t: np.array, m_t: np.array):

        exog = np.c_[s_t, np.repeat(1, len(s_t))]  # exog.shape => (nobs, 2)

        super(MultipleYsModel, self).__init__(
            endog=np.c_[i_t, m_t], exog=exog, k_states=2, initialization="diffuse"
        )

        self.ssm["design"] = np.zeros((self.k_endog, self.k_states, self.nobs))
        self.ssm["design", 0, 0, :] = s_t
        self.ssm["design", 1, 1, :] = 1

        # These have ok shape. Placeholders since I'm changing them
        # in the update() function
        self.ssm["selection"] = np.eye(self.k_states)
        self.ssm["transition"] = np.eye(self.k_states)

        # Dictionary of positions to names
        self.position_dict = OrderedDict(
            var_e1=1, var_e2=2, var_w1=3, var_w2=4, delta1=5, delta2=6, delta3=7
        )
        self.initial_values = starting_values
        self.positive_parameters = slice(0, 4)

    @property
    def param_names(self):
        return list(self.position_dict.keys())

    @property
    def start_params(self):
        """
        Initial values
        """
        # (optional) Use scale for var_e1 and var_e2 starting values
        params = np.r_[
            self.initial_values["var_e1"],
            self.initial_values["var_e2"],
            self.initial_values["var_w1"],
            self.initial_values["var_w2"],
            self.initial_values["delta1"],
            self.initial_values["delta2"],
            self.initial_values["delta3"],
        ]
        return params

    def transform_params(self, unconstrained):
        """
        If you need to restrict parameters
        For example, variances should be > 0
        Parameters maybe have to be within -1 and 1
        """
        constrained = unconstrained.copy()
        constrained[self.positive_parameters] = (
            constrained[self.positive_parameters] ** 2
        )
        return constrained

    def untransform_params(self, constrained):
        """
        Need to reverse what you did in transform_params()
        """
        unconstrained = constrained.copy()
        unconstrained[self.positive_parameters] = (
            unconstrained[self.positive_parameters] ** 0.5
        )
        return unconstrained

    def update(self, params, **kwargs):
        params = super(MultipleYsModel, self).update(params, **kwargs)

        # The following line is not needed (by default, this matrix is initialized by zeroes),
        # But I leave it here so the dimensions are clearer
        self["obs_intercept"] = np.repeat([np.array([0, 0])], self.nobs, axis=0).T

        self["obs_cov", 0, 0] = params[0]
        self["obs_cov", 1, 1] = params[1]

        self["state_cov"] = np.diag(params[2:4])

        # delta1, delta2, delta3
        self["transition", 0, 0] = params[4]
        self["transition", 0, 1] = params[5]
        self["transition", 1, 1] = params[6]

[ ]:

mod = MultipleYsModel(i_hat, s_t, m_hat)
res = mod.fit()

print(res.summary())

ボーナス: 高速なベイズ推定のための pymc3

このセクションでは、カスタムの状態空間モデルを pymc3 に簡単に組み込み、ベイズ法で推定する方法を示します。特に、この例では、No-U-Turn Sampler（NUTS）と呼ばれるハミルトニアンモンテカルロのバージョンを使った推定を示します。

基本的には、このノートブックに含まれるアイデアをコピーしているので、詳細についてはそちらを確認してください。

[ ]:

# Extra requirements
import pymc3 as pm
import theano
import theano.tensor as tt

モデルに暗黙的に含まれる尤度関数をtheanoに接続するために、いくつかのヘルパー関数を定義する必要があります。

[ ]:

class Loglike(tt.Op):

    itypes = [tt.dvector]  # expects a vector of parameter values when called
    otypes = [tt.dscalar]  # outputs a single scalar value (the log likelihood)

    def __init__(self, model):
        self.model = model
        self.score = Score(self.model)

    def perform(self, node, inputs, outputs):
        (theta,) = inputs  # contains the vector of parameters
        llf = self.model.loglike(theta)
        outputs[0][0] = np.array(llf)  # output the log-likelihood

    def grad(self, inputs, g):
        # the method that calculates the gradients - it actually returns the
        # vector-Jacobian product - g[0] is a vector of parameter values
        (theta,) = inputs  # our parameters
        out = [g[0] * self.score(theta)]
        return out


class Score(tt.Op):
    itypes = [tt.dvector]
    otypes = [tt.dvector]

    def __init__(self, model):
        self.model = model

    def perform(self, node, inputs, outputs):
        (theta,) = inputs
        outputs[0][0] = self.model.score(theta)

モデル1で使用したデータを再度シミュレートします。また、再度fitを行い、結果を保存して、得られるベイズ事後分布と比較します。

[ ]:

y_t, x_t, w_t, beta_x, beta_w = gen_data_for_model1()
plt.plot(y_t)

[ ]:

mod = TVRegression(y_t, x_t, w_t)
res_mle = mod.fit(disp=False)
print(res_mle.summary())

ベイズ推定

各パラメータの事前分布と、ドロー数、バーンインポイント数を定義する必要があります。

[ ]:

# Set sampling params
ndraws = 3000  #  3000 number of draws from the distribution
nburn = 600  # 600 number of "burn-in points" (which will be discarded)

[ ]:

# Construct an instance of the Theano wrapper defined above, which
# will allow PyMC3 to compute the likelihood and Jacobian in a way
# that it can make use of. Here we are using the same model instance
# created earlier for MLE analysis (we could also create a new model
# instance if we preferred)
loglike = Loglike(mod)

with pm.Model():
    # Priors
    intercept = pm.Uniform("intercept", 1, 10)
    var_e = pm.InverseGamma("var.e", 2.3, 0.5)
    var_x_coeff = pm.InverseGamma("var.x.coeff", 2.3, 0.1)
    var_w_coeff = pm.InverseGamma("var.w.coeff", 2.3, 0.1)

    # convert variables to tensor vectors
    theta = tt.as_tensor_variable([intercept, var_e, var_x_coeff, var_w_coeff])

    # use a DensityDist (use a lamdba function to "call" the Op)
    pm.DensityDist("likelihood", loglike, observed=theta)

    # Draw samples
    trace = pm.sample(
        ndraws,
        tune=nburn,
        return_inferencedata=True,
        cores=1,
        compute_convergence_checks=False,
    )

事後分布は、MLE推定とどのように比較されますか？

明らかにMLE推定値の周りにピークがあります。

[ ]:

results_dict = {
    "intercept": res_mle.params[0],
    "var.e": res_mle.params[1],
    "var.x.coeff": res_mle.params[2],
    "var.w.coeff": res_mle.params[3],
}
plt.tight_layout()
_ = pm.plot_trace(
    trace,
    lines=[(k, {}, [v]) for k, v in dict(results_dict).items()],
    combined=True,
    figsize=(12, 12),
)
最終更新日：2024年10月3日