狠狠撸

資料実装例 rstan コード
1
MCMC 実装例例 1 と例 2 のコード
#Lee-Carter モデルのパラメータにある値を与えてそれが真のパラメータであるとします。それから計算され
る死亡率から 2 項分布である死亡者数をサンプリングして、それが観測されたデータだとします。今回はそ
の観測されたデータを入力値として MCMC を行うことで真のパラメータをどれくらい精度よく推定できてい
るかを確認します。
#MCMC には stan を使って推定しています。Stan は R 上で rstan というライブラリを呼び出して使用してい
ます。以下の実行にはまず rstan をインストールする必要がありますが、cran から直接は無理でいくらか作
業が必要です。検索すれば出てくるので必要であれば別途確認してください。
################
##実装例 1 と実装例 2 （年度?年齢による固定効果のみのモデル）
##例 1 と例 2 では入力する観察データが違うだけでコード本体は同じです。
## 1950 年～2010 年の 50 歳～90 歳の Lee-Carter モデルの真のパラメータを設定し各年齢各年度の死亡デー
タを作成しました。データの生成は excel で行い、R 入力形式に変換して入力しました。
library(rstan)
runnum <- 103
####観察データの入力####
r_inp0 <- c(略)
r_inp1 <- c(略)
r_inp2 <- c(略)
r_inp3 <- c(略)
n_inp0 <- c(略)
n_inp1 <- c(略)
n_inp2 <- c(略)
n_inp3 <- c(略)
r_inp <- c(r_inp0, r_inp1, r_inp2, r_inp3)
n_inp <- c(n_inp0, n_inp1, n_inp2, n_inp3)
data <- list(

2
nxt = structure(n_inp,.Dim=c(41,61)),
rxt = structure(r_inp,.Dim=c(41,61))
)
####パラメータ初期値の入力####
ax_inp <- c(-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-
5,-5,-5,-5)
bx_inp <-
c(0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.
01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01)
kt_inp <-
c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0)
init <- list(
ax = structure(ax_inp,.Dim=c(41)),
bx = structure(bx_inp,.Dim=c(40)),
kt = structure(kt_inp,.Dim=c(60))
)
####stan コードの本体####
stancode <- '
data {
int rxt[41,61];
int nxt[41,61];
}
parameters {
vector[41] ax;
vector[40] bx;
vector[60] kt;

3
}
model {
real pxt;
real bxsum;
real ktsum;
bxsum <- 0;
ktsum <- 0;
ax ~ normal(0,1.0E3);
bx ~ normal(0,1.0E3);
bxsum <- sum(bx);
kt ~ normal(0,1.0E3);
ktsum <- sum(kt);
for (x in 1:40){
for (t in 1:60){
pxt <- fmin(1.0,exp(ax[x] + bx[x] * kt[t]));
rxt[x,t] ~ binomial(nxt[x,t],pxt);
}
}
for (t in 1:60){
pxt <- fmin(1, exp(ax[41] + (1-bxsum) * kt[t]));
rxt[41,t] ~ binomial(nxt[41,t],pxt);
}
pxt <- fmin(1, exp(ax[41] + (1-bxsum) * (0-ktsum)));
rxt[41,61] ~ binomial(nxt[41,61],pxt);
}
'

4
#コードのコンパイル
sc <- stanc(model_code = stancode, model_name="sc")
sm <- stan_model(stanc_ret=sc, varbose=FALSE)
####1000 本（うち最初の 500 本は捨て）の sampling####
fit <- sampling(sm, data=data, iter = 1000,
init=list(init), chains = 1, seed=1111)
result <- extract(fit)
assign(paste("result",runnum, sep=""),extract(fit))
#####ax のグラフ出力#####
par(new=F)
jpeg(paste("ax",runnum,".jpg"))
ymin=min(result$ax); ymax=max(result$ax);
for (i in 1:500){
plot(1:41,result$ax[i,],xlim=c(1,40),ylim=c(ymin,ymax),type="l",xlab="Age-50",ylab="ax")
par(new=T)
}
plot(1:41,seq(-5,-2,length=41),xlim=c(1,40),col="red",lwd=3,ylim=c(ymin,ymax),type="l",xlab="Age-
50",ylab="ax")
dev.off()
#####bx のグラフ出力#####
par(new=F)
jpeg(paste("bx",runnum,".jpg"))
ymin=min(result$bx); ymax=max(result$bx);
for (i in 1:500){
plot(1:40,result$bx[i,],xlim=c(1,40),ylim=c(ymin,ymax),type="l",xlab="Age-50",ylab="bx")
par(new=T)
}

5
plot(1:41,c(rep(0.035,8),rep(0.025,12),seq(0.025,0.015,length=21)),xlim=c(1,40),col="red",lwd=3,ylim=c(ymin,
ymax),type="l",xlab="Age-50",ylab="bx")
dev.off()
#####kt のグラフ出力#####
par(new=F)
jpeg(paste("kt",runnum,".jpg"))
ymin=min(result$kt); ymax=max(result$kt);
for (i in 1:500){
plot(1:60,result$kt[i,],xlim=c(1,60),ylim=c(ymin,ymax),type="l",xlab="Year-1950",ylab="kt")
par(new=T)
}
plot(1:61,seq(30,-30,length=61),xlim=c(1,60),col="red",lwd=3,ylim=c(ymin,ymax),type="l",xlab="Year-
1950",ylab="kt")
dev.off()

6
MCMC 実装例例 3 のコード
################
##実装例 3（年齢のみによる固定効果のほかに、UW の質の違いなど何らかの要因によるランダム効果を含
むモデル）
## 2008 年～2010 年の 50 歳～90 歳の Lee-Carter モデルの真のパラメータを設定し各年齢各年度の死亡デー
タを作成しました。データの生成は excel で行い、R 入力形式に変換して入力しました。
library(rstan)
runnum <- 8
####観察データの入力####
r_inp0 <- c(略)
r_inp1 <- c(略)
r_inp <- c(r_inp0, r_inp1)
data_ms <- list(
rxt = structure(r_inp,.Dim=c(41,3,10))
)
####パラメータの初期値####
ax_inp <- c(-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-5,-
5,-5,-5,-5,-5)
rdx_inp <- 1
rx_inp <- rep(0,10)
init_ms <- list(
ax = structure(ax_inp,.Dim=c(41)),
rdx = rdx_inp,
rx = structure(rx_inp,.Dim=c(10))
)

7
####stan コード本体####
stancode_mixed_s <- '
data {
int rxt[41,3,10];
//int nxt[41,3];
}
parameters {
vector[41] ax;
real <lower=0> rdx;
vector[10] rx;
}
model {
real pxt;
ax ~ normal(0,1.0E3);
rx ~ normal(0,rdx);
rdx ~ uniform(0,1.0E3);
for (x in 1:41){
for (t in 1:3){
for (i in 1:10){
pxt <- fmin(1.0,exp(ax[x]+rx[i]));
rxt[x,t,i] ~ binomial(100,pxt);
}
}
}

8
}
'
#コードのコンパイル
sc_ms <- stanc(model_code = stancode_mixed_s, model_name="sc_mixed")
sm_ms <- stan_model(stanc_ret=sc_ms, varbose=FALSE)
####1000 本（うち最初の 500 本は捨て）の sampling####
fit_ms <- sampling(sm_ms, data=data_ms, iter = 1000,
init=list(init_ms),chains = 1, seed=1111)
result <- extract(fit_ms)
assign(paste("result",runnum,sep=""),extract(fit_ms))
rxadj <- mean(result$rx)
#####ax のグラフ出力#####
par(new=F)
jpeg(paste("ax",runnum,".jpg"))
ymin=min(result$ax); ymax=max(result$ax);
for (i in 1:500){
plot(1:41,result$ax[i,]+mean(result$rx[i,]),xlim=c(1,40),ylim=c(ymin,ymax),type="l",xlab="Age-50",ylab="ax")
par(new=T)
}
#真のパラメータ
ax_true=c(-5.7,-5.625,-5.55,-5.475,-5.4,-5.325,-5.25,-5.175,-4.9,-4.825,-4.75,-4.675,-4.6,-4.525,-4.45,-4.375,-
4.3,-4.225,-4.15,-4.075,-4,-3.915,-3.83,-3.745,-3.66,-3.575,-3.49,-3.40499999999999,-3.32,-3.235,-
3.14999999999999,-3.06499999999999,-2.97999999999999,-2.89499999999999,-2.80999999999999,-
2.72499999999999,-2.63999999999999,-2.55499999999999,-2.46999999999999,-2.38499999999999,-
2.29999999999999)
plot(1:41,ax_true,xlim=c(1,40),col="red",lwd=3,ylim=c(ymin,ymax),type="l",xlab="Age-50",ylab="ax")
dev.off()

9
#####ri のグラフ出力#####
par(new=F)
jpeg(paste("rx",runnum,".jpg"))
ymin=min(result$rx-rxadj); ymax=max(result$rx-rxadj);
for (i in 1:500){
plot(1:10,result$rx[i,]-mean(result$rx[i,]),xlim=c(1,10),ylim=c(ymin,ymax),type="l",xlab="i",ylab="ri")
par(new=T)
}
#真のパラメータ
rx_true=c(0, -0.665040568, -0.392788633, -0.827113654, 1.268953693, -0.025035057,
1.151139571, -0.76772477, 0.656378306, -0.402994197)
plot(1:10,rx_true,xlim=c(1,10),ylim=c(ymin,ymax),type="l",col="red",lwd=3,xlab="i",ylab="ri")
dev.off()
#####rdx の出力#####
mean(result$rdx)
jpgname <- paste("hist_sigma_r.jpg")
jpeg(jpgname)
hist(result$rdx,nclass=50, main="Histgram of sigma_r",xlab="sigma_r")
dev.off()

10
リスク量計測部分のコード
###################
xn=36
tn=51
inu=1
dx1p<-rep(0,500)
dx1pp<-array(0,dim=c(500,10000))
dx1avg<-0
dx2p<-rep(0,500)
dx2avg<-0
dx3p<-rep(0,500)
dx3avg<-0
mx1<-rep(0,500)
mx1avg<-0
mx2<-rep(0,500)
mx2avg<-0
mx3<-rep(0,500)
mx3avg<-0
for (k in 1:500){
mx1[k] <- exp(result102$ax[k,xn]+result102$bx[k,xn]*result102$kt[k,tn])
mx2[k] <- exp(result103$ax[k,xn]+result103$bx[k,xn]*result103$kt[k,tn])
mx3[k] <- exp(result8$ax[k,xn]+result8$rx[k,inu])
}
mx1avg <- exp(mean(result102$ax[,xn]+result102$bx[,xn]*result102$kt[,tn]))
mx2avg <- exp(mean(result103$ax[,xn]+result103$bx[,xn]*result103$kt[,tn]))
mx3avg <- exp(mean(result8$ax[,xn]+result8$rx[,inu]))

11
dx1avg <- 10000*mx1avg
for (k in 1:500){
dx1p[k] <- 10000*mx1[k]
dx1pp[k,] <- rbinom(10000,10000, mx1[k])
dx2p[k] <- 100*mx2[k]
dx2pp[k,] <- rbinom(10000,100, mx2[k])
dx3p[k] <- 100*mx3[k]
dx3pp[k,] <- rbinom(10000,100, mx3[k])
}
res <- ""
res <- rbind(res,c(dx1avg/10000,quantile(dx1p[],.995)/10000,quantile(dx1pp[],.995)/10000
,dx2avg/100,quantile(dx2p[],.995)/100,quantile(dx2pp[],.995)/100
,dx3avg/100,quantile(dx3p[],.995)/100,quantile(dx3pp[],.995)/100
))
write(res, file="res2.csv") #res2.csv として結果セットが出力される

狠狠撸

「ベイズ推定でパラメータリスクを捉える＆優れたサンプラーとしてのMCMC」の実装例rstanコード

More Related Content

「ベイズ推定でパラメータリスクを捉える＆優れたサンプラーとしてのMCMC」の実装例rstanコード

狠狠撸

「ベイズ推定でパラメータリスクを捉える ＆優れたサンプラーとしてのMCMC」の実装例rstanコード

More Related Content

「ベイズ推定でパラメータリスクを捉える ＆優れたサンプラーとしてのMCMC」の実装例rstanコード

「ベイズ推定でパラメータリスクを捉える＆優れたサンプラーとしてのMCMC」の実装例rstanコード

「ベイズ推定でパラメータリスクを捉える＆優れたサンプラーとしてのMCMC」の実装例rstanコード