2015-03-09

Stanで混合正規分布の推定②

混合正規分布のパラメータ推定シリーズ。
これまでEMアルゴリズムとStanのMCMCにてパラメータの推定を行ってきた。

今回は、数値最適化によるパラメータ推定をやってみようと思う。
StanといえばHMCを用いたMCMCだが、実は最近のバージョンアップにより数値最適化関数が用意されたようなので、今回はそれを使ってみる。

前回に引き続き、データは以下の分布からサンプリングすることにより作る。

${ x \sim 0.7 \mathcal{N} (-5, 5^2) + 0.3 \mathcal{N} (5, 1^2) }$

ここに出てくる2つの平均、2つの分散、2つの混合比率を推定するのが目的。

サンプリングされたデータのヒストグラムは以下の通りで、乱数の種を揃えているので、いつもと完全に同じデータ。

f:id:lofas:20150305160039p:plain

数値最適化のやり方は非常に簡単で、Stanのモデルコードを読み込んだあとに、モデルのsamplingメソッドを使っていたところを、optimizingメソッドに変更するだけでOK。
Stanコードに書いたパラメータの制約条件もちゃんと考慮してくれる。

In [36]: fitchan = model.optimizing(data=stan_data)

Stanの数値最適化関数を実行することで得られた各パラメータの推定結果はこちら。
EMアルゴリズムとStanのMCMCの推定結果とほぼ一致していることが確認できる。

In [29]: fitchan    
Out[29]: OrderedDict([(u'mu', array([-5.23176606,  5.00483165])), (u'sig', array([ 4.82527209,  1.02056465])), (u'pi', array([ 0.69367853,  0.30632147]))])

ただし、Stanの最適化関数は特に何も指定しない限り初期値をランダムに選択するようで、初期値によっては全然違った結果が返ってくることもある。
以下は、想定の答えが返ってこなかった例。

In [24]: fitchan
Out[24]: OrderedDict([(u'mu', array([-2.09644408, -1.57536974])), (u'sig', array([ 6.22396311,  6.08568693])), (u'pi', array([  9.99260820e-01,   7.39179944e-04]))])

このように、初期値にかなり依存するので使う場合は注意が必要。
EMアルゴリズムも初期値に依存する点は同様だが、よく言われるように、その依存度合いは数値最適化よりは小さい印象。

Pythonコードはこちらで、sampling→optimizingの1行だけ変わっている。

# -*- coding: utf-8 -*-
from __future__ import print_function, division
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
import seaborn as sns
import pystan

#サンプル数
N = 5000
#混合数
K = 2
#混合係数
pi = 0.7
#乱数の種
np.random.seed(0)

#混合係数から各分布のサンプリング数を決める
N_k1 = N*pi
N_k2 = N-N_k1

#真のパラメータ
mu1 = -5
sig1 = np.sqrt(25)
mu2 = 5
sig2 = np.sqrt(1)

x1 = np.random.normal(mu1,sig1,N_k1)
x2 = np.random.normal(mu2,sig2,N_k2)

#観測変数
x = np.hstack((x1,x2))
base=np.linspace(np.min(x),np.max(x),1000)
plt.hist(x,bins=100,normed=True)
plt.plot(base,pi*stats.norm.pdf(base,mu1,sig1))
plt.plot(base,(1-pi)*stats.norm.pdf(base,mu2,sig2))
plt.plot(base,pi*stats.norm.pdf(base,mu1,sig1)+(1-pi)*stats.norm.pdf(base,mu2,sig2))

#Stan
stan_data = {'N': N, 'M': K, 'x': x}

model = pystan.StanModel('D:/Python/GMM.stan')

fitchan = model.optimizing(data=stan_data)

fitchan

Stanコードはこちらで、前回と全く同じ。

data {
  int<lower=1> N;
  int<lower=1> M;
  real x[N];
}

parameters {
  vector[M] mu;
  vector<lower=0.0001>[M] sig;
  simplex[M] pi;
}

model {
    real ps[M];
    for(n in 1:N){
        for(m in 1:M){
            ps[m] <- log(pi[m]) + normal_log(x[n], mu[m], sig[m]);
        }
        increment_log_prob(log_sum_exp(ps));
    }
}

2015-03-05

Stanで混合正規分布の推定①

前回↓に引き続き、混合正規分布のパラメータ推定。ただし、今回はEMアルゴリズムではなくStanを使ってMCMCにてパラメータを推定する。

<a href="http://lofas.hatenablog.com/entry/2015/03/03/153615">EMアルゴリズムで混合正規分布の推定 - WHAT' CHA GONNA DO FOR ME?</a>

前回に引き続き、データは以下の分布からサンプリングすることにより作る。

${ x \sim 0.7 \mathcal{N} (-5, 5^2) + 0.3 \mathcal{N} (5, 1^2) }$

ここに出てくる2つの平均、2つの分散、2つの混合比率を推定するのが目的。

サンプリングされたデータのヒストグラムは以下の通りで、乱数の種を揃えているので、前回と完全に同じデータとなる。
ただし今回は、密度関数の値はそれぞれの分布のものに加えて、それらの合計も紫の線で示した。

f:id:lofas:20150305160039p:plain

Stanを実行することで得られたパラメータ推定値の分布とパスのプロットがこちら。
上段がそれぞれの分布の平均、中段が標準偏差（分散ではない）、下段が混合係数となっていて、それぞれ問題なく分布収束しているように見える。
バーンインは100しか取ってないのにStanすごいやん。

f:id:lofas:20150305155256p:plain

続いて、パラメータ推定値のサマリー。
EMアルゴリズムの結果とほぼ完璧に一致していることが確認できる。
ただ、EMアルゴリズムはパラメータの標準誤差を計算するのが難しい（がんばればできるようだが）ことから信頼区間が構築しづらいのに対し、MCMCでは以下のように簡単に信用区間を構築できるため、その点はMCMCの方が便利。

In [18]: fitchan
Out[18]: 
Inference for Stan model: anon_model_912e3ad2e9344e95d4acd093ad7b9b1b.
1 chains, each with iter=2000; warmup=100; thin=1; 
post-warmup draws per chain=1900, total post-warmup draws=1900.

         mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat
mu[0]   -5.24  7.7e-3   0.11  -5.44  -5.31  -5.24  -5.17  -5.01  194.0    1.0
mu[1]     5.0  1.3e-3   0.03   4.94   4.98    5.0   5.02   5.06  511.0    1.0
sig[0]   4.82  3.2e-3   0.07   4.68   4.78   4.82   4.87   4.97  531.0    1.0
sig[1]   1.02  1.1e-3   0.03   0.97   1.01   1.02   1.04   1.07  488.0    1.0
pi[0]    0.69  4.1e-4 7.7e-3   0.68   0.69   0.69    0.7   0.71  365.0    1.0
pi[1]    0.31  4.1e-4 7.7e-3   0.29    0.3   0.31   0.31   0.32  365.0    1.0
lp__   -1.5e4    0.09   1.59 -1.5e4 -1.5e4 -1.5e4 -1.5e4 -1.5e4  341.0    1.0

Samples were drawn using NUTS(diag_e) at 03/05/15 15:46:37.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

Pythonコードはこちら。

# -*- coding: utf-8 -*-
from __future__ import print_function, division
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
import seaborn as sns
import pystan

#サンプル数
N = 5000
#混合数
K = 2
#混合係数
pi = 0.7
#乱数の種
np.random.seed(0)

#混合係数から各分布のサンプリング数を決める
N_k1 = N*pi
N_k2 = N-N_k1

#真のパラメータ
mu1 = -5
sig1 = np.sqrt(25)
mu2 = 5
sig2 = np.sqrt(1)

x1 = np.random.normal(mu1,sig1,N_k1)
x2 = np.random.normal(mu2,sig2,N_k2)

#観測変数
x = np.hstack((x1,x2))
base=np.linspace(np.min(x),np.max(x),1000)
plt.hist(x,bins=100,normed=True)
plt.plot(base,pi*stats.norm.pdf(base,mu1,sig1))
plt.plot(base,(1-pi)*stats.norm.pdf(base,mu2,sig2))
plt.plot(base,pi*stats.norm.pdf(base,mu1,sig1)+(1-pi)*stats.norm.pdf(base,mu2,sig2))

#Stan
stan_data = {'N': N, 'M': K, 'x': x}

model = pystan.StanModel('D:/Python/GMM.stan')

fitchan = model.sampling(data=stan_data, iter=2000, warmup=100,chains=1)

fitchan    

fitchan.plot()
plt.tight_layout()

Stanコードはこちら。

data {
  int<lower=1> N;
  int<lower=1> M;
  real x[N];
}

parameters {
  vector[M] mu;
  vector<lower=0.0001>[M] sig;
  simplex[M] pi;
}

model {
    real ps[M];
    for(n in 1:N){
        for(m in 1:M){
            ps[m] <- log(pi[m]) + normal_log(x[n], mu[m], sig[m]);
        }
        increment_log_prob(log_sum_exp(ps));
    }
}

2015-03-03

EMアルゴリズムで混合正規分布の推定

EMアルゴリズムの導入として、定番の混合正規分布のパラメータ推定をやってみる。

問題としては、K個の正規分布からある比率に従ってサンプリングされたデータだけが手元に得られた状態で、それぞれの正規分布の平均と分散、及びその比率を推定するというもの。

数式で書けば、混合正規分布は次のようになる。

${\displaystyle x \sim\ \sum_{k=1}^{K}\pi_k \mathcal{N} (\mu_k, \sigma_k^2), \ \sum_{k=1}^{K}\pi_k=1 }$

今回は、最も簡単な ${K=2}$ のケースを設定する。
具体的には、次の分布を使う。

${ x \sim 0.7 \mathcal{N} (-5, 5^2) + 0.3 \mathcal{N} (5, 1^2) }$

ここに出てくる2つの平均、2つの分散、2つの混合比率を推定するのが目的。

まずは、この分布からデータをサンプリングするとこんな感じになる。
理論的な密度関数の値をヒストグラムに重ねてみても大体それに沿っているので、多分想定通りにサンプリングできていると思われる。

f:id:lofas:20150303145322p:plain

このデータだけをもとに(=一つひとつのデータがそれぞれどちらの正規分布からサンプリングされたか知ることなしに)、EMアルゴリズムを使ってパラメータを推定しにいく。
EMアルゴリズムは、パラメータに適当な初期値を与えた上で、EステップとMステップを繰り返すことで（この繰り返す数をイテレーションと呼んだりする。そこはMCMCと同じ）、尤度の高いところのパラメータを見つけに行く最尤推定の一種で、この繰り返し計算により尤度は単調に増加することが知られている（ただし、大域的最適解に辿り着くとは限らない⇔局所解に陥る可能性あり）

EMアルゴリズムのアルゴリズムの説明をここで行うのはあまりに大変なのと、既存の本よりうまく説明できる気も全くしないので、興味がある方は例えば下記の本を参照頂きたい。いずれも混合正規分布を例にした説明がある。
今回はおもに一番上のPRMLを参考にしてコードは書いた。

パターン認識と機械学習下 (ベイズ理論による統計的予測)

作者: C.M.ビショップ,元田浩,栗田多喜夫,樋口知之,松本裕治,村田昇
出版社/メーカー: 丸善出版
発売日: 2012/02/29
メディア: 単行本
購入: 6人クリック: 14回
この商品を含むブログを見る

はじめてのパターン認識

作者: 平井有三
出版社/メーカー: 森北出版
発売日: 2012/07/31
メディア: 単行本（ソフトカバー）
購入: 1人クリック: 7回
この商品を含むブログ (2件) を見る

これなら分かる最適化数学―基礎原理から計算手法まで

作者: 金谷健一
出版社/メーカー: 共立出版
発売日: 2005/09/01
メディア: 単行本
購入: 29人クリック: 424回
この商品を含むブログ (41件) を見る

以下は、EMアルゴリズムを用いて実際に推定を行った結果。
気になるパラメータの前に、まずは対数尤度が本当に単調増加しているかを確認する。
大体イテレーションが100ぐらいのところまでは対数尤度は単調に増加し、それ以降は上昇はストップしている模様。

f:id:lofas:20150303145330p:plain

最後のイテレーションのパラメータを抜いてくると、このような値が得られる。
上段が平均、中段が分散、下段が混合比率を表している。微妙なズレはあるものの、ほぼ正解と一致していることがわかる。

In [21]: theta[:,:,iteration-1]
Out[21]: 
array([[ -5.23184557,   5.00480952],
       [ 23.28299505,   1.04155592],
       [  0.69366902,   0.30633098]])

コード全体はこちら。
本当であれば、対数尤度の増加がある閾値を下回ったところでアルゴリズムを打ち切りたいところだが、今回はイテレーション数を頭のところで決め打ちにしている。
また、今回はEMアルゴリズムでパラメータ推定を行ったが、MCMCを用いたパラメータ推定も可能と思われる。それについてはまた別の機会に書くかもしれないし書かないかもしれない。

# -*- coding: utf-8 -*-
from __future__ import print_function, division
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
import seaborn as sns

#サンプル数
N = 5000
#混合数
K = 2
#混合係数
pi = 0.7
#イテレーション
iteration = 500
#乱数の種
np.random.seed(0)

#混合具合を一様乱数から決める
#N_k1 = np.round(np.random.uniform()*N,decimals=0)
N_k1 = N*pi
N_k2 = N-N_k1

#真のパラメータ
mu1 = -5
sig1 = np.sqrt(25)
mu2 = 5
sig2 = np.sqrt(1)

x1 = np.random.normal(mu1,sig1,N_k1)
x2 = np.random.normal(mu2,sig2,N_k2)

#観測変数
x = np.hstack((x1,x2))
base=np.linspace(np.min(x),np.max(x),1000)
plt.hist(x,bins=100,normed=True)
plt.plot(base,pi*stats.norm.pdf(base,mu1,sig1))
plt.plot(base,(1-pi)*stats.norm.pdf(base,mu2,sig2))

#負担率格納用
gamma = np.zeros((N,K,iteration))
#尤度格納用
llk = np.zeros((iteration,1))
#パラメータ格納用
theta = np.zeros((3,K,iteration))
#初期値
theta[0,0,:]=0
theta[0,1,:]=1
theta[1,0,:]=1
theta[1,1,:]=2
theta[2,0,:]=0.6
theta[2,1,:]=0.4

for i in range(iteration):
  ##################################
  #E-Step
  ##################################
    gamma[:,0,i] = theta[2,0,i]*stats.norm.pdf(x,theta[0,0,i],np.sqrt(theta[1,0,i]))/ \
        (theta[2,0,i]*stats.norm.pdf(x,theta[0,0,i],np.sqrt(theta[1,0,i]))+ \
         theta[2,1,i]*stats.norm.pdf(x,theta[0,1,i],np.sqrt(theta[1,1,i])))
         
    gamma[:,1,i] = theta[2,1,i]*stats.norm.pdf(x,theta[0,1,i],np.sqrt(theta[1,1,i]))/ \
        (theta[2,0,i]*stats.norm.pdf(x,theta[0,0,i],np.sqrt(theta[1,0,i]))+ \
         theta[2,1,i]*stats.norm.pdf(x,theta[0,1,i],np.sqrt(theta[1,1,i])))    

  ##################################
  #M-Step
  ##################################
    if i != iteration-1:
        #mu
        theta[0,0,i+1] = gamma[:,0,i].dot(x)/np.sum(gamma[:,0,i])
        theta[0,1,i+1] = gamma[:,1,i].dot(x)/np.sum(gamma[:,1,i])
        #sig
        theta[1,0,i+1] = gamma[:,0,i].dot((x-theta[0,0,i+1])**2)/np.sum(gamma[:,0,i])
        theta[1,1,i+1] = gamma[:,1,i].dot((x-theta[0,1,i+1])**2)/np.sum(gamma[:,1,i])
        #pi
        theta[2,0,i+1] = np.sum(gamma[:,0,i])/N
        theta[2,1,i+1] = np.sum(gamma[:,1,i])/N
  ##################################
  #対数尤度の計算
  ##################################
        llk[i] = np.sum(np.log( \
              theta[2,0,i+1]*stats.norm.pdf(x,theta[0,0,i+1],np.sqrt(theta[1,0,i+1]))+ \
              theta[2,1,i+1]*stats.norm.pdf(x,theta[0,1,i+1],np.sqrt(theta[1,1,i+1]))  \
        ))  

plt.plot(llk[range(iteration-1)])
theta[:,:,iteration-1]

2015-02-27

PyStanでOLS①

今回は、流行りのStanをPythonから使ってみようと思う。PyStanを使うことでPythonからStanが呼べるようになる。
Stanは、ハミルトニアンモンテカルロ(HMC)とメトロポリス・ヘイスティングス(MH)の合わせ技で、従来のMHアルゴリズムによるMCMCより効率的なサンプリングが可能と言われている。
さらに、NUTS(No-U-Turn Sampler)という技が使われていて、これはHMCを実装する際にユーザーが決めるべきパラメータを自動で決めてくれる優れもの。
この辺は、理論的裏付けが下記の論文に書いてあるようなので、いつか真面目に勉強したい。
http://www.stat.columbia.edu/~gelman/research/published/nuts.pdf

Stanは、RやPython、Matlabから呼べて、RからStanを呼ぶRStanはそれなりに触ったことがあるが、PyStanは使ったことがないので、今回触ってみる。

推定に使うモデルはいつもの超簡単な単回帰式。

${ y=2+1.5x+\epsilon\\ x \sim\ \mathcal{N}(1,1)\\ \epsilon \sim\ \mathcal{N}(0,1)\\ }$

推定したいモデルは、Pythonコードとは別に用意したStanコードの方に書いておく。
Pythonコードからは、Stanコードにデータを渡してキックするイメージ。

PyStanを実行することで得られるパラメータ推定値の分布とパスのプロットがこちら。
上段が定数項と回帰係数、下段が残差の標準偏差となっている。定数項と回帰係数はベクトルにまとめたため二つ同時に出ている。

f:id:lofas:20150227174404p:plain

続いて、パラメータ推定値のサマリー。
95%信用区間に真のパラメータ（ベイズの考え方からは本来こういう言い方はしないが‥）が含まれていることがわかる。

In [105]: fitchan
Out[105]: 
Inference for Stan model: anon_model_b665a1feb6ebc5ee8386c19e8517f522.
1 chains, each with iter=2000; warmup=100; thin=1; 
post-warmup draws per chain=1900, total post-warmup draws=1900.

          mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat
beta[0]   2.04  2.3e-3   0.04   1.96   2.01   2.04   2.07   2.13  352.0    1.0
beta[1]   1.47  1.6e-3   0.03   1.41   1.45   1.47   1.49   1.53  381.0    1.0
sigma     0.97  9.5e-4   0.02   0.93   0.95   0.97   0.99   1.01  555.0    1.0
lp__    -468.7    0.07   1.25 -471.8 -469.2 -468.3 -467.8 -467.3  312.0    1.0

Samples were drawn using NUTS(diag_e) at 02/27/15 17:19:26.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

一応、OLSの結果とも比較しておくと、ほぼ完全に一致していることが確認できる。
Stanのパラメータの平均値と標準偏差は、OLSのパラメータ推定値と標準誤差にそれぞれ相当する。

print(results.summary())
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.692
Model:                            OLS   Adj. R-squared:                  0.692
Method:                 Least Squares   F-statistic:                     2241.
Date:                Fri, 27 Feb 2015   Prob (F-statistic):          2.28e-257
Time:                        17:27:16   Log-Likelihood:                -1386.1
No. Observations:                1000   AIC:                             2776.
Df Residuals:                     998   BIC:                             2786.
Df Model:                           1                                         
==============================================================================
                 coef    std err          t      P>|t|      [95.0% Conf. Int.]
------------------------------------------------------------------------------
const          2.0431      0.043     47.942      0.000         1.959     2.127
x1             1.4691      0.031     47.339      0.000         1.408     1.530
==============================================================================
Omnibus:                        1.076   Durbin-Watson:                   2.074
Prob(Omnibus):                  0.584   Jarque-Bera (JB):                1.155
Skew:                           0.059   Prob(JB):                        0.561
Kurtosis:                       2.883   Cond. No.                         2.53
==============================================================================

Pythonコードはこちら。

# -*- coding: utf-8 -*-
from __future__ import print_function, division
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns
import statsmodels.api as sm
import pystan

#サンプリング数
N=1000

#データ生成
np.random.seed(0)
x=np.random.normal(1,1,N)
eps=np.random.normal(0,1,N)
y=2+1.5*x+eps

#OLS(statsmodels)
X=sm.add_constant(x, prepend=True)
model = sm.OLS(y,X)
results = model.fit()
print(results.summary())

#Stan
N = X.shape[0]
K = X.shape[1]

stan_data = {'N': N, 'K': K, 'y': y, 'X': X}

model = pystan.StanModel('D:\Python/OLS.stan')
fitchan = model.sampling(data=stan_data, iter=2000, warmup=100,chains=1)
                  
fitchan.plot()
fitchan

Stanコードはこちら。

data {
  int<lower=1> N;
  int<lower=1> K;
  real y[N];
  matrix[N,K] X;
}

parameters {
  vector[K] beta;
  real<lower=0.0001> sigma;
}

model {
  y ~ normal(X*beta,sigma);
}

2015-02-25

カルマンスムージング(KS)

前回は、カルマンフィルターをやったので、今回はスムージング（平滑化）をやろうと思う。

<a href="http://lofas.hatenablog.com/entry/2015/02/21/144227">カルマンフィルタ - WHAT' CHA GONNA DO FOR ME?</a>

モデルは前回と同様、下記のシンプルなモデルとする。パラメータも同じ値をセットした。

${ \begin{eqnarray} y_t &=& x_t+\sigma_{\epsilon}\epsilon_t, \ \epsilon_t \sim\ \mathcal{N}(0,1), \nonumber \\ x_t &=& \phi x_{t-1}+\sigma_{\eta}\eta_t, \ \eta_t \sim\ \mathcal{N}(0,1). \nonumber \end{eqnarray} }$

前回は明示的に書かなかったが、フィルタリングにおいては、 ${t}$ 期の観測値 ${y_t}$ が得られた下での ${t}$ 期の潜在変数 ${x_t}$ の分布 ${p(x_t|y_t)}$ を求めにいった。

これに対して、スムージングにおいては、 ${T}$ 期までの観測値 ${y_{1:T}}$ がすべて得られたもとでの各 ${t}$ 期の潜在変数 ${x_t}$ の分布 ${p(x_t|y_{1:T})}$ が求められる。

フィルタリングに比べて、より多くの観測値の情報を用いて潜在変数を推定することになるため、一般的に精度はより高くなる。

アルゴリズムとしては、カルマンフィルターで ${1}$ 期から ${T}$ 期までの潜在変数の平均と分散を計算した後、それらの値を使いながら ${T}$ 期から ${1}$ 期まで戻ってくるようなイメージとなる。

for t in np.arange(T,0,-1)-1:        
    if t==T-1 :
        x_smtd[t] = x_fltd[t]
        P_smtd[t] = P_fltd[t]
    else:
        x_smtd[t] = x_fltd[t]+P_fltd[t]*phi/P_fcst[t+1]*(x_smtd[t+1]-x_fcst[t+1])
        P_smtd[t] = P_fltd[t]+(P_fltd[t]*phi/P_fcst[t+1])**2*(P_smtd[t+1]-P_fcst[t+1])

なぜこのようなアルゴリズムになるのか興味がある方は、前回挙げたHamiltonやTsayに加えて、以下の本も参考になると思います。
HamiltonやTsayはファイナンス色が強いので、あまり馴染みのない方はこちらの方がいいかもしれません。

データ同化入門 (予測と発見の科学)

作者: 樋口知之
出版社/メーカー: 朝倉書店
発売日: 2011/09/15
メディア: 単行本
購入: 2人クリック: 8回
この商品を含むブログ (1件) を見る

フィルタリングとスムージングそれぞれについて、平均と上下2標準偏差の信頼区間をプロットしておく。

まずは、フィルタリング。

f:id:lofas:20150225132521p:plain

次に、スムージング。

f:id:lofas:20150225132530p:plain

スムージングの方が、線がより滑らかになっていることが見て取れる。

コード全体はこちら。

# -*- coding: utf-8 -*-
from __future__ import print_function, division
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd

T=200
phi=0.7
sig_eps=1
sig_eta=1

np.random.seed(0)
x=np.zeros((T,1))
y=np.zeros((T,1))

#潜在変数の作成
for t in range(T-1):
    x[t+1]=phi*x[t]+np.random.normal(0,sig_eta,1)
#観測値の作成
for t in range(T):
    y[t]=x[t]+np.random.normal(0,sig_eps,1)
    
#フィルタリングに必要な変数の箱
x_fcst=np.zeros((T,1))
P_fcst=np.zeros((T,1))
y_fcst=np.zeros((T,1))
MSE_y=np.zeros((T,1))
x_fltd=np.zeros((T,1))
P_fltd=np.zeros((T,1))
x_smtd=np.zeros((T,1))
P_smtd=np.zeros((T,1))

#潜在変数の平均と分散の初期値
x_fcst[0] = 0
P_fcst[0] = 1000

#カルマンフィルター
for t in range(T):
    
    y_fcst[t] = x_fcst[t]
    MSE_y[t] = P_fcst[t]+sig_eps**2
    
    x_fltd[t] = x_fcst[t]+P_fcst[t]/MSE_y[t]*(y[t]-y_fcst[t])
    P_fltd[t] = P_fcst[t]-P_fcst[t]**2/MSE_y[t]
    
    if t != T-1 :
        x_fcst[t+1] = phi*x_fltd[t]
        P_fcst[t+1] = phi**2*P_fltd[t]+sig_eta**2

#カルマンスムージング        
for t in np.arange(T,0,-1)-1:        
    if t==T-1 :
        x_smtd[t] = x_fltd[t]
        P_smtd[t] = P_fltd[t]
    else:
        x_smtd[t] = x_fltd[t]+P_fltd[t]*phi/P_fcst[t+1]*(x_smtd[t+1]-x_fcst[t+1])
        P_smtd[t] = P_fltd[t]+(P_fltd[t]*phi/P_fcst[t+1])**2*(P_smtd[t+1]-P_fcst[t+1])

#Kalman Filter
plt.plot(x_fltd-2*np.sqrt(P_fltd),lw=1,color="b")
plt.plot(x_fltd,lw=1,color="k")
plt.plot(x_fltd+2*np.sqrt(P_fltd),lw=1,color="b")
plt.title("Kalman Filter",fontsize=16)

#Kalman Smoothing
plt.plot(x_smtd-2*np.sqrt(P_smtd),lw=1,color="b")
plt.plot(x_smtd,lw=1,color="k")
plt.plot(x_smtd+2*np.sqrt(P_smtd),lw=1,color="b")
plt.title("Kalman Smoothing",fontsize=16)

2015-02-21

カルマンフィルタ(KF)

ある時系列データが観測値として得られているとして、これらの時系列データがどのようにして生成されているか、その仕組みに興味があるとする。
ただ、いま手元にあるのは観測値として得られた時系列データだけなので、その背後にある仕組みについては何かしらの方法で推定する必要がある。
状態空間モデルというフレームワークを用いることで、この推定が可能となる。

今回は、次のようなシンプルな状態空間モデルを考えてみたい。

ここで、 ${y_t}$ は ${t}$ 期の観測値、 ${x_t}$ は ${t}$ 期の観測値を駆動する潜在変数とする。
${x_t}$ がAR(1)モデルに従いながら背後で動いていて、この挙動が興味の対象であるものの、我々が観測できるのは ${x_t}$ にノイズが加わった ${y_t}$ だけ、ということをこのモデルは言っている。

上記の式のうち、上の式は観測方程式、下の式は状態方程式と呼ばれており、これらのセットが状態空間モデルである。
観測方程式、状態方程式が共に線形で、かつノイズが正規分布に従う場合に限り、カルマンフィルタというアルゴリズムにより潜在変数 ${x_t}$ の推定は効率的に行われることが知られている。

というわけで、前置きが長くなったが、今回はシミュレーションデータを使ってカルマンフィルタをPythonで書いてみたよ、というお話。

まず、観測値 ${y_t}$ 、潜在変数 ${x_t}$ のそれぞれを以下のように作る。

In [19]: for t in range(T-1):
    ...:     x[t+1]=phi*x[t]+np.random.normal(0,sig_eta,1)
    ...: 
    ...: #観測値の作成
    ...: for t in range(T):
    ...:     y[t]=x[t]+np.random.normal(0,sig_eps,1)
    ...: 
    ...: 
    ...: #観測値と潜在変数のプロット
    ...: plt.plot(np.hstack((y,x)))
    ...: #plt.title("Observed Value vs Latent Variable")
    ...: plt.legend(['Observed Value','Latent Variable'])

Out[19]: <matplotlib.legend.Legend at 0x11d0f8f0>

f:id:lofas:20150221142701p:plain

潜在変数に対してノイズが加わっている観測値は、その分変動が大きいことが見て取れる。

なお、モデルの未知パラメータである ${\phi}$ 、 ${\sigma_{\epsilon}}$ 、 ${\sigma_{\eta}}$ は、それぞれ0.7、1、1をセットした。

次はメインのカルマンフィルタ部分。

In [29]: for t in range(T):
    ...:     
    ...:     y_fcst[t] = x_fcst[t]
    ...:     MSE_y[t] = P_fcst[t]+sig_eps**2
    ...:     
    ...:     x_fltd[t] = x_fcst[t]+P_fcst[t]/MSE_y[t]*(y[t]-y_fcst[t])
    ...:     P_fltd[t] = P_fcst[t]-P_fcst[t]**2/MSE_y[t]
    ...:     
    ...:     if t != T-1 :
    ...:         x_fcst[t+1] = phi*x_fltd[t]
    ...:         P_fcst[t+1] = phi**2*P_fltd[t]+sig_eta**2

なぜこのようなアルゴリズムになるのか興味がある方は、例えば以下の本を参照してください。

Time Series Analysis

作者: James D. Hamilton
出版社/メーカー: Princeton Univ Pr
発売日: 1994/01/11
メディア: ハードカバー
購入: 1人クリック: 5回
この商品を含むブログ (8件) を見る

Analysis of Financial Time Series (Wiley Series in Probability and Statistics)

作者: Ruey S. Tsay
出版社/メーカー: Wiley
発売日: 2010/08/30
メディア: ハードカバー
クリック: 5回
この商品を含むブログを見る

推定された潜在変数と真の潜在変数とを比べてみる。

f:id:lofas:20150221142723p:plain

まあ大体トラックできているような感じはする。
この差分も念のためプロットしておく。

f:id:lofas:20150221142739p:plain

最後に、観測値とも一緒にプロットしてみる。

f:id:lofas:20150221142755p:plain

観測値に含まれるノイズに騙されることなく、よりマイルドな動きをする真の潜在変数に近い値が推定できていることが、なんとなく見て取れる。

今回はデータをシミュレーションから作ったため、パラメータはあらかじめこちらで与えた。
しかし、実際の応用では観測データしか手元に無いため、潜在変数とあわせてパラメータも推定の対象となることが多い。
これについては、また別の機会に書くかもしれないし書かないかもしれない。

コード全体はこちら。

# -*- coding: utf-8 -*-
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd

T=200
phi=0.7
sig_eps=1
sig_eta=1

np.random.seed(0)
x=np.zeros((T,1))
y=np.zeros((T,1))

#潜在変数の作成
for t in range(T-1):
    x[t+1]=phi*x[t]+np.random.normal(0,sig_eta,1)
#観測値の作成
for t in range(T):
    y[t]=x[t]+np.random.normal(0,sig_eps,1)
    
#観測値と潜在変数のプロット
plt.plot(np.hstack((y,x)))
#plt.title("Observed Value vs Latent Variable")
plt.legend(['Observed Value','Latent Variable'])

#フィルタリングに必要な変数の箱
x_fcst=np.zeros((T,1))
P_fcst=np.zeros((T,1))
y_fcst=np.zeros((T,1))
MSE_y=np.zeros((T,1))
x_fltd=np.zeros((T,1))
P_fltd=np.zeros((T,1))

#潜在変数の平均と分散の初期値
x_fcst[0] = 0
P_fcst[0] = 1000

#カルマンフィルター
for t in range(T):
    
    y_fcst[t] = x_fcst[t]
    MSE_y[t] = P_fcst[t]+sig_eps**2
    
    x_fltd[t] = x_fcst[t]+P_fcst[t]/MSE_y[t]*(y[t]-y_fcst[t])
    P_fltd[t] = P_fcst[t]-P_fcst[t]**2/MSE_y[t]
    
    if t != T-1 :
        x_fcst[t+1] = phi*x_fltd[t]
        P_fcst[t+1] = phi**2*P_fltd[t]+sig_eta**2

#真の潜在変数と推定された潜在変数のプロット
plt.plot(np.hstack((x,x_fltd)))
plt.legend(['True Latent','Estimated Latent'])

#推定された潜在変数の真の潜在変数との差分
plt.plot(x-x_fltd)

#真の潜在変数と推定された潜在変数とプロット
plt.plot(np.hstack((x,x_fltd,y)))
plt.legend(['True Latent','Estimated Latent','Observed Value'])

2015-02-20

pandasで日次リターンを月次リターンに変換する②

日次リターンを月次リターンに変換する技、第2弾。

今回もこちらのありがたい本で紹介されている方法を踏襲する。

Pythonによるデータ分析入門 ―NumPy、pandasを使ったデータ処理

作者: Wes McKinney,小林儀匡,鈴木宏尚,瀬戸山雅人,滝口開資,野上大介
出版社/メーカー: オライリージャパン
発売日: 2013/12/26
メディア: 大型本
この商品を含むブログ (7件) を見る

以下の2つの方法が紹介されているが、

①日次の累積リターンから計算する方法
②日次のリターンから計算する方法

今回は②の方法で、日次リターンを月次リターンに変換する。

いま、stock_retという変数に日次のリターンがすでに入っているとする。

stock_ret.plot()
Out[12]: <matplotlib.axes._subplots.AxesSubplot at 0x150de850>

f:id:lofas:20150219161511p:plain

これを使って日次リターンを月次リターンに変換するが、コードは次のようにするだけでOK。

ret_m=(1+stock_ret).resample("M",how="prod",kind="period")-1
ret_m.plot()
Out[13]: <matplotlib.axes._subplots.AxesSubplot at 0x15309b70>

f:id:lofas:20150219161432p:plain

Mは、month frequencyのことで、これでデータの頻度を日次から月次に変換する。
prodは、サマリー対象のレコード全てを掛け算した値を使う、という意味で、ここまでで各月の月次のグロスリターン系列ができることとなる。
最後に1を引くことで、ネットの月次リターンとなる。
timestampかperiodかを選べるが、今回はperiodを選択。

コード全体はこちら。

# -*- coding: utf-8 -*-
import pandas.io.data as web
import numpy as np
import seaborn as sns
import pandas as pd
import pandas.io.data as web
import datetime

start = datetime.datetime(2000, 1, 1)
end = datetime.datetime(2014, 12, 31)

#データ取得
all_data = {}
for ticker in ['AAPL', 'IBM', 'MSFT', 'YHOO']:
    all_data[ticker] = web.DataReader(ticker,'yahoo',start,end) 

#データ結合(concatでもOK)
#stock_price=all_data['AAPL'].ix[:,['Adj Close']]
#for ticker in ['IBM', 'MSFT', 'YHOO']:
#    stock_price = pd.concat((stock_price,all_data[ticker].ix[:,['Adj Close']]),axis=1)

#データ結合
stock_price=all_data['AAPL'].ix[:,['Adj Close']]
for ticker in ['IBM', 'MSFT', 'YHOO']:
    stock_price = pd.merge(stock_price,all_data[ticker].ix[:,['Adj Close']],
                           left_index=True, right_index=True)
#列名変更
stock_price.columns=['AAPL', 'IBM', 'MSFT', 'YHOO']
stock_price.plot()

#リターンデータに変換
stock_ret=stock_price.pct_change()
stock_ret.plot()

#月次リターンに変換
ret_m=(1+stock_ret).resample("M",how="prod",kind="period")-1
ret_m.plot()