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

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

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

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

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

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

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

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

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

<pre>
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
</pre>

<pre>
import matplotlib.pyplot as plt

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

![ダウンロード (15).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/163445/c6046510-08be-31cd-2370-5e769638000d.png)

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

<pre>
import pandas as pd
import numpy as np
import hmmlearn
from hmmlearn.hmm import GaussianHMM
</pre>

<pre>
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)
</pre>

![ダウンロード (17).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/163445/8b142391-5e99-691c-1a0b-59671ee374aa.png)

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

<pre>
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)
</pre>

<pre>
hidden_states, mus, sigmas, transmat, model = fitHMM(y, 2)
</pre>

<pre>
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()
</pre>

<pre>
t = np.array(range(len(y)))
</pre>

<pre>
plot_states(y, hidden_states, t)
</pre>

![ダウンロード (18).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/163445/ccecb981-fa4a-bda7-59b8-012d99f84885.png)

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

<pre>
hidden_states, mus, sigmas, transmat, model = fitHMM(y, 3)
</pre>

<pre>
plot_states(y, hidden_states, t)
</pre>

![ダウンロード (19).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/163445/662f5779-05c9-3a2c-d7b0-786f76cd5d6d.png)


https://github.com/RuoAndo/time-series-analysis-with-Python/blob/master/qiita/boston-hmm-1.ipynb