渋谷駅前で働くデータサイエンティストのブログ

元祖「六本木で働くデータサイエンティスト」です / 道玄坂→銀座→東京→六本木→渋谷駅前

Stanで統計モデリングを学ぶ(3): ざっと「Stanで何ができるか」を眺めてみる

R BUGS/Stan 統計学 時系列分析

実は業務でもStan使い始めてるんですが、まだまだ単位根ありパネルデータの分析に回すなど低レベルなものが多く、無情報事前分布と階層事前分布を巧みに使いこなして華麗にサンプリング。。。なんて夢のまた夢という情けない状況です（泣）。


で、気が付いたら@さんのStan関連ブログ記事が超絶充実していて、久保先生もびっくりみたいな状況に。もはや僕が何かをだらだら書くのもアホらしいので、先にStanの使い方を覚えたいという方は是非@さんのブログから読んで下さい（笑）。僕はひたすらそちらの記事を（例えばinfer.NETあたりの例題を解きながら）トレースしていくだけのショボい記事をだらだら書いていこうと思ってます。


ということで、とりあえずStanマニュアルであるstan-reference-2.1.0.pdfを斜め読みして*1、ざっと僕が興味のある範囲でまとめただけのものを書き出してみます。Stanコード例は全てマニュアルからの引用です。


はじめに


基本的にStan（というかMC / MCMCサンプラー）は尤度計算をモンテカルロ法に基づいてやってくれるものです。なので、尤度計算さえ収束すれば良いというタイプの統計モデリングは事実上何でもStanでできるはずです。よって、ここに名前が挙がってないモデリング手法も、うまくコードを書いて実装できればStanで走らせることが可能だとも言えると思います。


f:id:TJO:20140207180929p:plain


モンテカルロ法による尤度計算のイメージはこんな感じ、というのは前回記事で紹介した通りです。で、大事なことはどんなモデリングにせよ「元のモデル式をきちんと『左辺＝右辺（確率分布）』の形に直せる」ということ。例えば正規線形モデルであれば、


Y_n=\alpha+\beta x_n + \eps_n \hspace{10pt} where \hspace{10pt} \eps_n \sim Normal(0, \sigma)
↓
Y_n - (\alpha+\beta x_n) \sim Normal(0, \sigma)
↓
Y_n \sim Normal(\alpha+\beta x_n, \sigma)
↓

y[i] ~ normal( alpha + beta * x[i], sigma )


というようにStanコードで表現できるところまで式変形できるということですね。慣れればパッと書けるようになると思いますが、複雑なモデリングを行いたい場合は最初にベタッと原理に従ったモデリング式を書き、それを式変形してStanコードに改めるというのが良さそうです。


多変量解析まわり


大体何でもできますが、Stanはコード文法がややこしかったり、サンプリングの仕方によって収束の度合いがまるっきり変わったりするので注意が必要です。分からない時は迷わず@さんあたりに聞きましょうｗ（何でも丸投げｗｗ）


正規線形モデル


まず超の字がつく基本中の基本として、正規線形モデル。これは正規分布normal(mu, sigma)で計算できます。

＜例＞

data {
	int<lower=0> N;
	vector[N] x;
	vector[N] y;
}
parameters {
	real alpha;
	real beta;
	real<lower=0> sigma;
}
model {
	for (n in 1:N)
		y[n] ~ normal(alpha + beta * x[n], sigma);
}

ロジスティック回帰・プロビットモデル


これはベルヌーイ分布bernoulli(theta)で計算できます。目的変数が二値でちゃんとintで与えておかないと普通にコンパイルエラーになるので要注意。

＜例＞
data {
	int<lower=0> N;
	real x[N];
	int<lower=0,upper=1> y[N];
}
parameters {
	real alpha;
	real beta;
}
model {
	for (n in 1:N)
		y[n] ~ bernoulli(inv_logit(alpha + beta * x[n]));
}

多項ロジットモデル


二値ロジットモデルができるなら、当然多項ロジットもできます。これは多項分布categorical(theta)で実現可能です。ちなみに順序ロジスティック回帰はそのまんまのordered_logistic(eta,c)でできます。

＜例＞
data {
	int K;
	int N;
	int D;
	int y[N];
	vector[D] x[N];
}
parameters {
	matrix[K,D] beta;
}
model {
	for (k in 1:K)
		for (d in 1:D)
			beta[k,d] ~ normal(0,5);
	for (n in 1:N)
		y[n] ~ categorical(softmax(beta * x[n]));
}

その他


階層ロジスティックモデルとかitem-response theory model(1PL-Rasch model / Multi-level 2PL model)とかも出てくるんですが、僕はよく分からないのでここでは割愛します。。。*2


計量時系列分析まわり


何でも勝手に確率分布をくっつけてモデリングできるというStanの特性を生かして、普通に計量時系列モデルのパラメータ推定を行うことができます。「Rで計量時系列分析」シリーズ記事で見てきたように、これまた最尤法で推定しているのでもちろんStanでバッチリやれるというわけです。


ARモデル


タイムラグの分だけ差分を取って、正規分布normal(mu,sigma)でモデリングすればOKです。なお発展としてARCHモデルも推定できますが、GARCHモデルの説明が詳しいのでここでは割愛。

＜例：AR(1)モデル＞
data {
	int<lower=0> N;
	real y[N];
}
parameters {
	real alpha;
	real beta;
	real sigma;
}
model {
	for (n in 2:N)
		y[n] ~ normal(alpha + beta*y[n-1], sigma);
}

＜例：AR(K)モデル＞
data {
	int<lower=0> K;
	int<lower=0> N;
	real y[N];
}
parameters {
	real alpha;
	real beta[K];
	real sigma;
}
model {
	for (n in (K+1):N) {
		real mu;
		mu <- alpha;
		for (k in 1:K)
			mu <- mu + beta[k] * y[n-k];
		y[n] ~ normal(mu, sigma);
	}
}

GARCHモデル


実はARCH(1)モデルの説明に使われているパラメータ群の説明があるので要注意ｗ

＜例：GARCH(1,1)モデル＞
data {
	int<lower=0> T;
	real r[T];
	real<lower=0> sigma1;
}
parameters {
	real mu;
	real<lower=0> alpha0;
	real<lower=0,upper=1> alpha1;
	real<lower=0,upper=(1-alpha1)> beta1;
}
transformed parameters {
	real<lower=0> sigma[T];
	sigma[1] <- sigma1;
	for (t in 2:T)
		sigma[t] <- sqrt(alpha0
				+ alpha1 * pow(r[t-1] - mu, 2)
				+ beta1 * pow(sigma[t-1], 2));
}
model {
	r ~ normal(mu,sigma);
}

MAモデル


MA(Q)モデルはベクトル表現を使って以下のように書けます。

＜例：MA(Q)モデル＞
data {
	int<lower=0> Q; // num previous noise terms
	int<lower=3> T; // num observations
	vector[T] y; // observation at time t
}
parameters {
	real mu; // mean
	real<lower=0> sigma; // error scale
	vector[Q] theta; // error coeff, lag -t
}
transformed parameters {
	vector[T] epsilon; // error term at time t
	for (t in 1:T) {
		epsilon[t] <- y[t] - mu;
		for (q in 1:min(t-1,Q))
			epsilon[t] <- epsilon[t] - theta[q] * epsilon[t - q];
	}
}
model {
	vector[T] eta;
	mu ~ cauchy(0,2.5);
	theta ~ cauchy(0,2.5);
	sigma ~ cauchy(0,2.5);
	for (t in 1:T) {
		eta[t] <- mu;
		for (q in 1:min(t-1,Q))
			eta[t] <- eta[t] + theta[q] * epsilon[t - q];
	}
	y ~ normal(eta,sigma);
}

ARMAモデル


ここではARMA(1,1)モデルの例だけ引用しておきます。

＜例：ARMA(1,1)モデル＞
data {
	int<lower=1> T; // num observations
	real y[T]; // observed outputs
}
parameters {
	real mu; // mean coeff
	real phi; // autoregression coeff
	real theta; // moving avg coeff
	real<lower=0> sigma; // noise scale
}
model {
	vector[T] nu; // prediction for time t
	vector[T] err; // error for time t
	nu[1] <- mu + phi * mu; // assume err[0] == 0
	err[1] <- y[1] - nu[1];
	for (t in 2:T) {
		nu[t] <- mu + phi * y[t-1] + theta * err[t-1];
		err[t] <- y[t] - nu[t];
	}
	mu ~ normal(0,10); // priors
	phi ~ normal(0,2);
	theta ~ normal(0,2);
	sigma ~ cauchy(0,5);
	err ~ normal(0,sigma); // likelihood
}

その他


確率的ボラティリティモデルや隠れマルコフモデルをコードできるそうですが、後者はともかく前者は僕は普段は全然扱わなくて疎いので割愛しました。。。


測定誤差とメタアナリシス


この辺僕はまーーーーーーったく詳しくないので全て割愛（ごめんなさい）。とは言え、測定誤差のモデリングもStanを使えば自由自在にできるので、例にも挙がっているようにメタアナリシスなどで測定誤差同士を統合して評価したいケースなんかでは物凄く有用な気がします。


クラスタリングまわり


わざわざ多項分布categorical(theta)が実装されているのを見れば分かる通り、ベイジアンの流儀でクラスタリングを実装することもできます。


K-means


ごくごく普通のユークリッド距離を使ったバージョンがマニュアルに載っています。

＜例："Soft" K-means＞
data {
	int<lower=0> N; // number of data points
	int<lower=1> D; // number of dimensions
	int<lower=1> K; // number of clusters
	vector[D] y[N]; // observations
}
transformed data {
	real<upper=0> neg_log_K;
	neg_log_K <- -log(K);
}
parameters {
	vector[D] mu[K]; // cluster means
}
transformed parameters {
	real<upper=0> soft_z[N,K]; // log unnormalized clusters
	for (n in 1:N)
		for (k in 1:K)
		soft_z[n,k] <- neg_log_K
				- 0.5 * dot_self(mu[k] - y[n]);
}
model {
	// prior
	for (k in 1:K)
		mu[k] ~ normal(0,1);
	// likelihood
	for (n in 1:N)
	increment_log_prob(log_sum_exp(soft_z[n]));
}

「ベイジアンで生成モデルやるのは難しいんだよコラ」


と、何故か1節まるまる使ってマニュアルで物凄く説教してます（笑）。詳しくはマニュアル本文を読んでもらいたいのですが、Non-IdentifiabilityとMultimodalityが元凶であるというコメントがされてます。ひとまずここは@さんにでもご解説を賜るとして*3、次にいきましょう。


混合モデル


実はマニュアルのトップに出てる例なんですが、普通に混合正規分布とかも多項分布categorial(theta)と正規分布normal(mu,sigma)との合わせ技でいけます。

＜例：混合モデル＞
data {
	int<lower=1> K; // number of mixture components
	int<lower=1> N; // number of data points
	real y[N]; // observations
}
parameters {
	simplex[K] theta; // mixing proportions
	real mu[K]; // locations of mixture components
	real<lower=0,upper=10> sigma[K]; // scales of mixture components
}
model {
	real ps[K]; // temp for log component densities
	for (k in 1:K) {
		mu[k] ~ normal(0,10);
	}
	for (n in 1:N) {
		for (k in 1:K) {
			ps[k] <- log(theta[k])
				+ normal_log(y[n],mu[k],sigma[k]);
		}
		increment_log_prob(log_sum_exp(ps));
	}
}

Latent Dirichlet Allocation (LDA)


ということで色々制約はあるらしいんですが、事前分布にディリクレ分布dirichlet(alpha)を充てて、尤度計算を多項分布categorical(theta)でやれば、ズバリLatent Dirichlet Allocationを走らせることができます。

＜例：LDA＞
data {
	int<lower=2> K; // num topics
	int<lower=2> V; // num words
	int<lower=1> M; // num docs
	int<lower=1> N; // total word instances
	int<lower=1,upper=V> w[N]; // word n
	int<lower=1,upper=M> doc[N]; // doc ID for word n
	vector<lower=0>[K] alpha; // topic prior
	vector<lower=0>[V] beta; // word prior
}
parameters {
	simplex[K] theta[M]; // topic dist for doc m
	simplex[V] phi[K]; // word dist for topic k
}
model {
	for (m in 1:M)
		theta[m] ~ dirichlet(alpha); // prior
	for (k in 1:K)
		phi[k] ~ dirichlet(beta); // prior
	for (n in 1:N) {
		real gamma[K];
		for (k in 1:K)
			gamma[k] <- log(theta[doc[n],k]) + log(phi[k,w[n]]);
		increment_log_prob(log_sum_exp(gamma)); // likelihood
	}
}

その他


ガウス過程を使うやつとか、cholesky_decompose()関数でコレスキー分解が出来るとか*4、色々なトピックスが他にも続くんですが、僕の現在の理解の範疇を今回は一旦割愛します。また必要になったら取り上げますよーということで。むしろこれからのStan修行大変かも。。。

*1:まだStan2.2.0にアップデートしてないんですごめんなさい

*2:勉強しろよゴルァとか言わないでー

*3:マジこればっかｗ

*4:これ知らなかった。。。