Scikit-learn Boston Housing Datasetを隠れマルコフモデルで解析してみる

■　マルコフ性とは
日本の気温の変化は、ヨーロッパの気温の変化に影響を及ぼしにくい。
→　システム（地球）全体ではなく、近傍を見るだけでよい、というのがマルコフ性の考え方。

■　隠れマルコフ性とは
天気を知りたいのだが、気温しか観測できない場合。

晴れた日は気温が高くなる可能性が高い
雨の日は気温が低くなる可能性が高い

しかし、測定できるのは気温だけ。

このような場合、隠れマルコフモデルを使う。

例えば、「晴れ」の隠れ状態は気温が高くなる観測を生み出すことがあるが、なんらかの理由で気温が低くなる観測がされることがある。

①曇り（気温：普通）→②雨（気温：普通）→③雨（気温：寒い）の場合、２の雨の隠れ状態では気温が普通である可能性もある。

前置きが長くなったが、ボストンハウジングデータセットを使ってみる。。

from sklearn import preprocessing
import pandas as pd
import numpy as np # we'll need it later
# Load the Boston dataset.

from sklearn.datasets import load_boston

boston = load_boston()
X,y = boston.data, boston.target

import matplotlib.pyplot as plt

plt.figure(figsize=(24, 8), dpi=50)
plt.subplot(211)
plt.plot(X)
plt.subplot(212)
plt.plot(y)

今回は、ターゲットデータのｙだけを使う。

import pandas as pd
import numpy as np
import hmmlearn
from hmmlearn.hmm import GaussianHMM

vals = np.expand_dims(y, 1)
n_states = 2
model = GaussianHMM(n_components=n_states, n_iter=100).fit(vals)
hidden_states = model.predict(vals)

下の部分が隠れマルコフの状態で、例えば、景気が良い／悪いなどの状態かも知れない。

def fitHMM(vals, n_states):
    vals = np.reshape(vals,[len(vals),1])
    
    # fit Gaussian HMM to Q
    model = GaussianHMM(n_components=n_states, n_iter=100).fit(vals)
     
    # classify each observation as state 0 or 1
    hidden_states = model.predict(vals)
 
    # fit HMM parameters
    mus = np.squeeze(model.means_)
    sigmas = np.squeeze(np.sqrt(model.covars_))
    transmat = np.array(model.transmat_)
    print(mus)
    print(sigmas)
    
#     relabeled_states = [state_dict[h] for h in hidden_states]
    relabeled_states = hidden_states
    return (relabeled_states, mus, sigmas, transmat, model)

hidden_states, mus, sigmas, transmat, model = fitHMM(y, 2)

def plot_states(ts_vals, states, time_vals):
    
    fig, ax1 = plt.subplots(figsize=(24.0, 8.0))
    color = 'tab:red'
    ax1.set_xlabel('time')
    ax1.set_ylabel('traffic', color=color)
    ax1.plot(time_vals, ts_vals, color=color)
    ax1.tick_params(axis='y', labelcolor=color)

    ax2 = ax1.twinx()  
    color = 'tab:blue'
    ax2.set_ylabel('Hidden state', color=color)  
    ax2.plot(time_vals,states,     color=color)
    ax2.tick_params(axis='y', labelcolor=color)

    fig.tight_layout()  
    plt.figure(figsize=(24, 8), dpi=50)
    plt.show()

t = np.array(range(len(y)))

plot_states(y, hidden_states, t)

さらに状態数を３つにしてみる。。。

hidden_states, mus, sigmas, transmat, model = fitHMM(y, 3)

plot_states(y, hidden_states, t)