rstan
- ようやくrstanがわかってきたのでまとめ直す
- まとめ
- データがある
- 確率モデルを定め、対数尤度計算をコード化する
- rのc++連携の仕組みを使ったrstanパッケージのstan()関数を使ってStan推定をする
- 何がよいのか
- 上記のまとめは読んでもつまらないまとめになっている
- だが、何も書いていないことが、Stanよさ
- 対数尤度計算ができさえすれば、データと確率モデルに制約はない
- この線で、単純なところからやってみる
- データが無い場合
- 推定という意味では意味がない設定ではあるが、事前分布設定をしているのだ、ということを確認するという意味でやってみる
- 何かしらの現象が確率pで起きているものとする。事前確率として一様分布を採用する
- pという変数名を設定する。それはreal型であって、下限値0、上限値1で定める
- model{}の中身は：
- 対数尤度をloglikという名前で計算することにする。初めに値0を代入する。0-1の一様分布から、あるpが得られる確率を考える。それはある観察pの0-1一様分布の尤度であって、それは、pが0-1の間にある限り一定であるが、それの対数(対数尤度)をloglikに加算している。increment_log_prob()はStanにおいて大事な関数で、最終的な対数尤度を、Stan的に記録、利用する処理を行っている
- nodata.stanというファイル名で以下を作る
functions{}
data{}
transformed data{}
parameters{
real<lower=0,upper=1> p;
}
transformed parameters{}
model{
real loglik;
loglik <- 0;
loglik <- loglik + uniform_log(p,0,1);
increment_log_prob(loglik);
}
generated quantities{}
-
- 実行すると、事前分布として指定した一様分布が推定結果として返る
out1 <- stan(file="nodata.stan") stan_hist(out1)
- データが１個の場合
- stanファイルのdata{}にデータを受け取れるように書き加える
- たった１個の0か1の値を与えるだけなので、下限値を0とした整数変数を宣言する
functions{}
data{
int<lower=0> x;
}
transformed data{}
parameters{
real<lower=0,upper=1> p;
}
transformed parameters{}
model{
real loglik;
loglik <- 0;
loglik <- loglik + uniform_log(p,0,1);
increment_log_prob(loglik);
}
generated quantities{}
-
- 実行に当たっては
x <- list(x=1) out2 <- stan(file="onedata.stan",data=x) stan_hist(out2)
- データが計画的に得られている場合
- N回観察してn回起きたとする
- データとして、長さNの0,1のベクトルがあるとすれば、要素数Nとベクトルをdataとして渡し、要素ごとに対数尤度を加算していけばよい
functions{}
data{
int<lower=0> n;
int<lower=0> x[n];
}
transformed data{}
parameters{
real<lower=0,upper=1> p;
}
transformed parameters{}
model{
real loglik;
loglik <- 0;
loglik <- loglik + uniform_log(p,0,1);
for(i in 1:n){
if(x[i]==0){
loglik <- loglik + bernoulli_log(0,p);
}else{
loglik <- loglik + bernoulli_log(1,p);
}
}
increment_log_prob(loglik);
}
generated quantities{}
-
- 以下でも同じこと
functions{}
data{
int<lower=0> n;
int<lower=0> x[n];
}
transformed data{}
parameters{
real<lower=0,upper=1> p;
}
transformed parameters{}
model{
real loglik;
loglik <- 0;
loglik <- loglik + uniform_log(p,0,1);
for(i in 1:n){
loglik <- loglik + bernoulli_log(x[i],p);
}
increment_log_prob(loglik);
}
generated quantities{}dat <- list(x=c(0,0,1,1,1,0,1),n=length(x)) out3 <- stan(file="ndata.stan",data=dat) stan_hist(out3)
-
- loglikを計算するのはいつものことなので、特定の分布から標本が得られた尤度を計算してincrement_log_probを加算していく過程を、"~"を使って次のように書きなおす方法がある
functions{}
data{
int<lower=0> n;
int<lower=0> x[n];
}
transformed data{
}
parameters{
real<lower=0,upper=1> p;
}
transformed parameters{}
model{
p ~ uniform(0,1);
for(i in 1:n){
x[i] ~ bernoulli(p);
}
}
generated quantities{}dat <- list(x=c(0,0,1,1,1,0,1),n=length(x)) out4 <- stan(file="ndata2.stan",data=dat)
-
- 対数尤度の計算でループを回さずに、ベクトル演算で一気にやるには
functions{}
data{
int<lower=0> n;
int<lower=0> x[n];
}
transformed data{
}
parameters{
real<lower=0,upper=1> p;
}
transformed parameters{}
model{
p ~ uniform(0,1);
x ~ bernoulli(p);
}
generated quantities{}out5 <- stan(file="ndata3.stan",data=dat)
- データが雑多な場合、その１
- m事象ずつまとめて観察した。まとめた観察をn回行ったところ、n回中k回は、mのうちすべてで成功し失敗は0だったという
- 成功確率pのとき、m個中m個とも成功する確率は
なので、そのようなモデルで書いてみる
functions{}
data{
int<lower=0> n;
int<lower=0> m;
int<lower=0> k;
}
transformed data{
}
parameters{
real<lower=0,upper=1> p;
}
transformed parameters{
real q;
q <- p^m;
}
model{
p ~ uniform(0,1);
k ~ binomial(n,q);
}
generated quantities{}
-
- 実行すると、推定対象はpだが、pから派生させてパラメタqを定義し、modelはqを使って書いている
- 推定はpとqとで行ってある
- pとqとは１対１対応だが、その分布の形は異なることに注意
out6 <- stan(file="irreg1.stan",data=list(n=10,m=3,k=4)) stan_hist(out6)
- データが雑多な場合、その２
- m事象ずつまとめて観察した。まとめた観察をn回行ったところ、n回中k回は、mのうちすべてで成功し失敗は0だったという
- ただし、n回のまとめ観察の内訳は、ある時はm1事象ずつまとめ、ある時はm2事象ずつまとめてあり、記録内容は、mi個中、全mi個が成功だったかどうかという0/1の情報だという
functions{}
data{
int<lower=0> n;
int<lower=0> m[n];
int<lower=0> k[n];
}
transformed data{
}
parameters{
real<lower=0,upper=1> p;
}
transformed parameters{
}
model{
p ~ uniform(0,1);
for(i in 1:n){
k[i] ~ bernoulli(p^m[i]);
}
}
generated quantities{}
-
- 実行
m <- c(2,5,3,2,3,3,3,6,7,3) k <- c(0,0,1,1,0,0,1,1,0,0) out7 <- stan(file="irreg2.stan",data=list(n=length(m),m=m,k=k))
- データが雑多な場合、その３
- もっといい加減な記録だったとする
- m事象ずつまとめて観察した。まとめた観察をn回行ったところ、n回中k回は、mのうちすべてで成功し失敗は0だったという
- ただし、n回のどれでmi中miだったかの情報はなく、単に、mという長さnの整数列と、kという整数しかないようなときに、どうしたらよいだろう
- k/nという観察は、１観察なので、尤度を計算するのは１回きりになるだろう
- n回中k回、何かが起きたことがわかっている。その内訳は、
,p^m[2],...,p^m[n]]の確率で起きることが、起きたり起きなかったりしながら、n個のうちk個で起きているので、1,2,...,nのうちk個取り出してそれで起きていて、残りで起きていないなら、p^m[i]と(1-p^m[j])とを組み合わせについて掛け合わせた確率をすべてのk個の組み合わせについて計算する必要がある
- 以下のようにもできるが…
- これだったら、Rでqを出してから、それがn回中k回起きたときのqの尤度を計算してもよい…のでつまらないが、まあ、思考訓練と思ってやってみる
functions{}
data{
int<lower=0> n;
int<lower=0> m[n];
int<lower=0> k;
int a;
int comb[k,a];
}
transformed data{
}
parameters{
real<lower=0,upper=1> p;
}
transformed parameters{
real q;
real pm[n];
real negpm[n];
real prodneg;
real tmp;
for(i in 1:n){
pm[i] <- p^m[i];
negpm[i] <- 1-pm[i];
}
prodneg <- prod(negpm);
q <- 0;
for(i in 1:n){
tmp <- prodneg;
for(j in 1:k){
tmp <- tmp * pm[comb[j,i]]/negpm[comb[j,i]];
}
q <- q + tmp;
}
}
model{
p ~ uniform(0,1);
k ~ binomial(n,q);
}
generated quantities{}
- データが雑多な場合、その４
- m事象ずつまとめて観察した。まとめた観察をn回行ったところ、n回中k回は、mのうちすべてで成功し失敗は0だったという
- この事象はまとめて観察するときは、そのm事象について一定の生起確率があると考えられるが
- 別のタイミングで観察するときには、生起確率が異なってくるものとする
- そして、そのタイミングによって異なる生起確率は、どこかに平均があるがばらついているわけで、それがベータ分布になっているものと仮定して良いとする
- ベータ分布のパラメタ２つは1から大きな値までの一様分布を事前分布とすることとする
- ベータ分布のパラメタを推定するとともに、そこから計算できるタイミングでばらつく生起確率の平均値と分散との事後分布を作ってみる
functions{}
data{
int<lower=0> n;
int<lower=0> m[n];
int<lower=0,upper=1> k[n];
}
transformed data{
}
parameters{
real<lower=1> a;
real<lower=1> b;
real<lower=0,upper=1> p[n];
}
transformed parameters{
}
model{
a ~ uniform(1,100000);
b ~ uniform(1,100000);
p ~ beta(a,b);
for(i in 1:n){
k[i] ~ bernoulli(p[i]^m[i]);
}
}
generated quantities{
real ave;
real v;
ave <- (a+1)/(a+b+2);
v <- a*b/((a+b)^2*(a+b+1));
}m <- c(2,5,3,2,3,3,3,6,7,3) k <- c(0,0,1,1,1,0,1,0,0,1) n <- length(m) out9 <- stan(file="irreg4.stan",data=list(n=n,m=m,k=k)) stan_hist(out9,pars=c("a","b","ave","v")) stan_scat(out9,pars=c("a","b")) stan_scat(out9,pars=c("ave","v"))
- やりたいことをやってみる
- 複合ジェノタイプ、標本混合がStanに乗るかどうかを考えてみる
