상태 공간 모형을 설치해 봤어요.

제조업 출신 데이터 과학자가 보낸 글

이번에는 상태 공간 모형을 설치해 보았다.

개시하다

며칠 전에 소개한 시간 시퀀스 분석를 찾아보고 상태 공간 모델의 기법을 알고 실시해 봤습니다.

상태 공간 모델

지난번에 소개한 타임라인 모델은 현지에서 예처리하기 불편하면 사용할 수 없어 현실적으로 사용하기 편하지 않나.이럴 때 사용하는 수법은 바로 상태 공간 모델이다.
이번에는 상세한 이론을 정리한 것이 아니기 때문에 상태 공간 모델의 맛을 아래에 기재한다.

비선형 행위를 처리할 수 있는 시간 서열

관찰치가 정적 분포를 제외한 분포 상황을 처리할 수 있다

상태 공간 모형의 생각

상태 공간 모델의 기본적인 생각은 다음과 같다.

어느 시점에서의 관찰값은 해당 시점의 실제 상태에서 오차가 발생한다

특정 시점의 실제 상태는 한 시점 이전의 실제 상태와 유사하다

즉, 실제 상태는 시간에 따라 무작위로 변화한다

관측 방정식: 관측 값의 시간 서열

상태 방정식: 실제 상태의 시간 서열

python의 코드는 다음과 같습니다.

# 必要なライブラリーのインポート
import numpy as np
import pandas as pd
from numpy.random import *
from scipy import stats

import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline

# グラフを横長にする
from matplotlib import rcParams
rcParams['figure.figsize'] = 10, 6
sns.set()

# pystanの読み込み
import pystan

# データの読み込み
# https://stat.ethz.ch/R-manual/R-devel/library/datasets/html/AirPassengers.html
df = pd.read_csv('AirPassengers.csv')

# float型に変換
df['#Passengers'] = df['#Passengers'].astype('float64')
df = df.rename(columns={'#Passengers': 'Passengers'})

# datetime型にしてインデックスにする
df.Month = pd.to_datetime(df.Month)
df = df.set_index("Month")

# データの中身を確認
df.head()

# データの可視化
plt.plot(df.earnings)

이번에 사용된 데이터는 매달 비행기 승객 수 데이터다.대상 기간은 1949년 1월부터 1960년 12월까지다.
처음 나타난 것은 "Box, G.E.P., Jenkins, G.M.and Reinsel, G.C.(1976) Time Series Analysis, Forecting and Control.Theird Edition.Hold-Day. Series G."
URL에서 데이터의 상세한 내용을 확인할 수 있다.
그런 다음 Stan에서 로컬 레벨 모델을 추측합니다.

\begin{align*}
y_t &= \mu_t + \epsilon_t, &\epsilon_t \sim N(0, \sigma_y^2) \\
\mu_t &= \mu_{t-1} + \eta_t,  &\eta_t \sim N(0, \sigma_{\mu}^2)
\end{align*}

stan 코드는 다음과 같다.

local_level = '''
data {
    int<lower=0> T; // number of learning points
    int<lower=0> M; // number of predict points
    real Y[T]; // observations
}

parameters {
    real mu[T]; // trend
    real<lower=0> s_y; // sd of observations
    real<lower=0> s_mu; // sd of trend
}

transformed parameters {
    real y_hat[T]; // prediction

    for(t in 1:T) {
        y_hat[t] = mu[t];
    }
}

model {
    for(t in 1:T) {
        Y[t] ~ normal(mu[t], s_y);
    }
    for(t in 2:T) {
        mu[t] ~ normal(mu[t-1], s_mu);
    }
}

generated quantities {
    real mu_pred[T+M];
    real y_pred[T+M];

    for(t in 1:T) {
        mu_pred[t] = mu[t];
        y_pred[t] = y_hat[t];
    }
    for(t in (T+1):(T+M)) {
        mu_pred[t] = normal_rng(mu_pred[t-1], s_mu);
        y_pred[t] = mu_pred[t];
    }
}
'''

이어서 stan 파일의 컴파일을 진행합니다.

stan_model = pystan.StanModel(model_code=local_level)

그런 다음 다양한 설정을 수행합니다.

# 目的変数の指定
y = df["Passengers"]

# 学習期間と予測期間の時点を指定
T = 130 #学習期間
M = 14 #予測期間

# 目的変数を学習用と予測用に分割
y_train = y[:-M]
y_test = y[-M:]

다음은 stan에게 데이터를 건네주십시오.이때 데이터는 분류기를 통해 대상을 대입하는 것을 주의하십시오.

predict_dat = {'T': T, 'M' : M, 'Y': y_train}

그리고 Sampling 함수를 사용하여 뒷분포부터 샘플링을 시작합니다.

fit_local_level = stan_model.sampling(data=predict_dat, iter=3000, chains=1, seed=10, n_jobs=1)

샘플링 결과를 추출하다.Rhat은 후분포 수렴 여부를 판단하는 지표다.1.1 이하면 수렴으로 판단할 수 있다.

fit_local_level

그리고 샘플링 결과를 추출하여 결과를 그립니다.

# サンプリング結果の抽出
ms_local_level = fit_local_level.extract()

# 予測値
y_pred = ms_local_level['y_pred'].mean(axis=0)

quantile = [5, 95]
per_5_95 = np.percentile(ms_local_level['y_pred'], q=quantile, axis=0).T
colname = ['p5', 'p95']
df_pred = pd.DataFrame(per_5_95, columns=colname)

# 予測値を追加
df_pred['y_pred'] = y_pred
mu_hat = ms_local_level['mu'].mean(axis=0)

# 状態の推定値を追加
df_pred['mu_hat'] = np.nan
df_pred.loc[0:129,'mu_hat'] = mu_hat

df.plot(y="Passengers", legend=False) # 目的変数
plt.plot(df_pred[['p5','p95']][-14:], linestyle="dashed", color='purple') # 予測区間
plt.plot(df_pred[['y_pred']][-14:], color='red') # 予測値
plt.plot(mu_hat, color='green') # 状態
plt.show()

전망치는 대체로 직선을 그리며 신뢰 구간이 점차 확대되고 있다.
이 모형이 안 좋아서 이렇게 된 거야.
사실 여기서부터 개량 모델을 통해 적합한 모델을 구축하면 예측이 잘 되지만 다음 기사에서 정리하고 싶습니다.

최후

끝까지 읽어주셔서 감사합니다.
이번에는 pstan을 사용하여 상태 공간 모델을 실현해 보았다.
앞으로는 모델 개량 등이 필요하지만, 실제 데이터에 맞춰 실시됐으면 좋겠다고 생각한다.
수정 요청이 있으면 연락 주세요.