読者です 読者をやめる 読者になる 読者になる

廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

Stan の integrate_ode_rk45 を使ってバスモデルのパラメータ推定

R Stan 微分方程式

以下の微分方程式で記述されるモデルのパラメータを推定します。

 \displaystyle dx_{t}/dt=\left(a+bx_{t}\right)(1-x_t)

ほんとうは閉じた形で解が求まるのですが今回は Stan の integrate_ode_rk45 を使って数値的に解を求めます。

カラーテレビの普及率のデータ(第1章第2節3 1 情報通信機器の世帯普及率 : 平成18年版 情報通信白書)に対して当てはめます。

最尤推定

モデルを記述した Stan コードは以下です。正規分布を仮定してます。

//Bass.stan
functions {
  real[] Bass(real t,        // time
            real[] y,      // state
            real[] theta,  // parameters
            real[] x_r,    // data (real)
            int[] x_i) {   // data (integer)
    real dydt[1];
    dydt[1] = (theta[1]+theta[2]*y[1])*(1-y[1]);
    return dydt;
  }
}
data{
  int<lower=0> T;
  real<lower=0> TS[T];
  real<lower=0> y[T];
}
transformed data {
   real x_r[0];
   int x_i[0];
   real Y0[1];
   int T0;
   Y0[1] = 0;
   T0 = 0;
}
parameters{
  real<lower=0> theta[2];
  real<lower=0> sigma;
}
transformed parameters {
   real y_hat[T,1];
   y_hat = integrate_ode_rk45(Bass, Y0, T0, TS, theta, x_r, x_i);
}
model{
  y ~ normal(y_hat[,1],sigma);
}

R から Stan を呼びます。

library(rstan)
library(cowplot)
y <- c(0.3, 1.6, 5.4, 13.9, 26.3, 42.3, 61.1, 75.8, 85.9, 90.3,
       93.7, 95.4, 97.7, 97.8, 98.2, 98.5, 98.9, 98.8, 99.2, 99.1,
       98.9, 98.7, 99.0, 99.3, 99.4, 99.3, 99.0, 99.1, 99.0, 98.9,
       99.1, 99.2, 99.2, 98.9, 99.0, 99.2, 99.3, 99.4, 99.0, 99.3, 99.4)/100
year=1965:2005
dat <- data.frame(year,y)
n=length(y)
dat4stan <- list(y=y,T=n,TS=1:n)
rstan_options(auto_write = TRUE)
options(mc.cores = 4)
Bass <- stan_model("~/Documents/Bass.stan")
optBass <- optimizing(Bass,data=dat4stan)
ggplot(dat,aes(x=year,y=y))+
  geom_point()+
  geom_line(aes(x=year,y=optBass$par[-c(1:3)]),colour="royalblue",size=1)

f:id:abrahamcow:20161118050636p:plain

ベイズ推定

パラメータには無情報の一様事前分布を仮定します。

最尤推定で求めた値を初期値に使います。

R のコードはこちら。

ini <- list(theta=optBass$par[1:2],sigma=optBass$par[3])
inits <-lapply(1:4, function(i)ini)
sampBass <- sampling(Bass,data=dat4stan,init=inits)
traceplot(sampBass,pars=c("theta","sigma"))
yhat <-get_posterior_mean(sampBass,pars="y_hat")[,"mean-all chains"]
ex <- rstan::extract(sampBass)
upper <-apply(ex$y_hat, 2, quantile,0.98)
lower <-apply(ex$y_hat, 2, quantile,0.01)
ggplot(dat,aes(x=year,y=y))+
  geom_point()+
  geom_ribbon(aes(ymin=lower,ymax=upper))+
  geom_line(aes(x=year,y=yhat),colour="royalblue",size=1)

前になにかでためしたときはめちゃめちゃ時間かかった記憶があるのですが、今回は一分未満で計算が終わりました。

f:id:abrahamcow:20161118050958p:plain

黒い帯で99%信頼区間を合わせてプロットしてあるのですが青い線(事後期待値)に重なっていてほとんど見えません。

トレースプロットはこちらです。

f:id:abrahamcow:20161118053317p:plain

ばっちり収束しています。

参考にしたい文献

微分方程式の解放には以下のアルゴリズムを使っているとのこと。

  • rk45: a fourth and fifth order Runge-Kutta method for non-stiff systems (Dormand and Prince, 1980; Ahnert and Mulansky, 2011)
    • Dormand, J. R. and Prince, P. J. (1980). A family of embedded Runge-Kutta formulae. Journal of Computational and Applied Mathematics, 6(1):19–26. 201
    • Ahnert, K. and Mulansky, M. (2011). Odeint—solving ordinary differential equations in C++. arXiv, 1110.3397.
  • bdf: a variable-step, variable-order, backward-differentiation formula implementation for stiff systems (Cohen and Hindmarsh, 1996; Serban and Hind- marsh, 2005)
    • Cohen, S. D. and Hindmarsh, A. C. (1996). CVODE, a stiff/nonstiff ODE solver in C. Computers in Physics, 10(2):138–143.
    • Serban, R. and Hindmarsh, A. C. (2005). CVODES: the sensitivity-enabled ODE solver in SUNDIALS. In ASME 2005 International Design Engineering Technical Confer- ences and Computers and Information in Engineering Conference, pages 257–269. American Society of Mechanical Engineers.

abrahamcow.hatenablog.com