StataでMCMC

前から興味はあったのですが、専門が専門なんで殆どベイズ統計学を学ばずに今まで過ごして来ました。しかし、MCMCの隆盛を見るにつけ、どうしても学びたくなってしまいました。大森先生のHPから色々と教材を頂いて来て、ざっくり初歩の初歩を学習。もっと敷居が高いと思っていたけど、今のところ思ったより簡単。これなら簡単なモデルならプログラムかけるぞと思い、Oxで書かれたサンプルプログラム（normal.ox）を、使い慣れたStataに移植してみました。計算してみると...Stata遅い...です。これは使えないかな。Oxで1秒かかってないのに、21.78秒もかかってしまいました。走らせている環境が違うので単純には比較できないとは思いますが、遅すぎます。MCMCはおとなしくOxでやっていた方が無難かな。

*追記（6/17）今日同じプログラムをまわしてみたら、7.1秒でした!? Stataの名誉回復の為に、もう少し補足しておきますと、Stata10はWindowsバージョンで、Macのエミュレーター（VMFusion）上で走っています。なので多少(?)ハンディ有りという環境です。昨日はそのエミュレーター上のWinXPの挙動がおかしかったからな。

// MCMC_NORMAL.do for STATA BY TOMBOYA on June 16, 2009
clear

timer clear 
timer on 1
// True parameters
scalar MUT=5
scalar VART=1

local n=100 // # obs
local r=1e4 // # repetitions
local nburn = 1e3 // # burn-ins

set obs `=`r''

// Data ganeration
gen x = rnormal(MUT,sqrt(VART)) in 1/`n'
qui sum x 
scalar XBAR = r(mean)
scalar SUMSQDEVX = r(Var) * (`n'-1)

// prameters of prior
// mu|var ~ N(mu0,var0)
// var ~ InvGamma(a0,b0)

// prior parameter for mu
scalar MU0=0
scalar VAR0=1000

// prior parameter for var
scalar A0=1
scalar B0=1

// variables storing outcomes
gen double sqdevx = . // square of deviation of x from mu
gen double mu_ = .
gen double var_ = .

// initial guess of true parameters
scalar MU = XBAR
scalar VAR = SUMSQDEVX/(`n'-1)

// draw (mu,var) from conditional posterior function
forvalue i = -`nburn'(1)`r'{
  // generate mu(i) given var(i-1) and x
  scalar VAR1 = 1/ (1/VAR0 + `n'/VAR)
  scalar MU1 = VAR1*(MU0/VAR0 + `n'*XBAR/VAR)
  scalar MU = rnormal(MU1,sqrt(VAR1))

  // generate var(i) given mu(i) and x using InvGamma(a,b)
  // mean=a/b, variance=a/b^2
  qui replace sqdevx = (x -MU)^2 in 1/`n'
  qui sum sqdevx
  scalar SUMSQDEVX = `n' * r(mean) 
  scalar A = A0 +`n'
  scalar B = B0 + 0.5*SUMSQDEVX
  // note that if x ~ Gamma(A,1/B), then 1/x ~ InvGamma(A,B) 
  scalar VAR = 1 / rgamma(A,1/B)

  if `i'>0{
    if mod(`i',1000)==0  dis "t=`i'"
    qui replace mu_ = MU in `i'
    qui replace var_= VAR in `i'
  }
}
timer off 1
timer list 1