トップ 差分 一覧 ソース 検索 ヘルプ PDF RSS ログイン

Bayes

Bayes

ベイズ法に関連するプログラム等。

最終更新時間：2009年12月11日 13時23分28秒


ベイズ法について


ベイズ法または、伝統的な統計学と異なり統計的推論にベイズの定理を使う統計学の手法です。古くから、優れた手法であることは知られていましたが、推定に必要な積分の計算が困難だったことから、統計学の主流とはなりませんでした。しかし、近年、コンピュータの処理速度の進歩とMCMC法の確立により、多くのモデルが推定可能となり，さまざまな分野への応用が進んでいます。

パッケージソフト


 BUGS


ベイズ法で事後分布を積分する際に必要となるモンテカルロ積分を実行するソフトです。MCMC法のうちギブス・サンプラーを使っています。いまのこころ、もっとも多様なモデルに対応しているようです。

BUGS系のソフトのうち、Windows用のものがWinBUGSです。


  • BUGSコードについては下記参照

 RのMCMCpackパッケージ


BUGSを使わなくても、Rの内部である程度のことはできます。他にも、bayesmパッケージやMNPパッケージがあります。


 パッケージソフトは使わない


もともと、パッケージソフトを使う目的は、時間の短縮のためです。統計学の方法の理解のためには、自分でプログラムを組むという作業が不可欠です。また、マルコフ切り替えモデルなど複雑なモデルの推定を行うには、パッケージソフトが対応していないので、自分でプログラミングを行う必要があります。

プログラミングには、一般的に MATLAB/Gauss が使われているようです。ただし、これらはインタープリタ型なので、スピードは遅めです。

Rの組み込み関数だけでMCMCを適用すると、スピードはかなり遅くなります。C/C++ と組み合わせて使う必要があります。


BUGSコード


 実例：株価の変化率の t 分布への当てはめ




データセットと一緒にワーキングディレクトリにおいておきます。

model{
for(i in 1:n){
e[i]~dt(mu,tau,nu)
}
mu ~dnorm(0, tau2)
tau ~dgamma(5,4)
tau2 <- 10*tau
nu ~dgamma(2, 0.5) 
}

  • R プログラム

#ライブラリ"R2WinBUG"の読み込み。事前にインストールしておく必要があります。
library(R2WinBUGS)
#データセットの読み込み。
topix <- read.csv("topix.csv",header=F)

topix.ts<-ts(topix[[2]])
e<-c(100*diff(log(topix.ts)))
n<-NROW(e)
in1<-list(mu=0,tau=2,nu=10)
in2<-list(mu=0,tau=2,nu=15)
in3<-list(mu=0,tau=2,nu=7.5)
inits<-list(in1,in2,in3)
data<-list("n","e")
parameters <- c("mu", "tau", "nu")

#MCMCの実行
topix.sim<-bugs(data,inits=inits, parameters, model.file="model_stock.txt",
n.chains = 3, n.iter = 5000, debug=F,
bugs.directory = "c:/Program Files/WinBUGS14/",
working.directory = NULL)

#結果の表示
print(topix.sim,digits=3)

#グラフ
mu <- topix.sim$mean$mu
tau <- topix.sim$mean$tau
df <- topix.sim$mean$nu
temp <- c(-10000:10000)/1000

hist(e,breaks=50,probability=1)
lines(temp/tau+mu,dt(temp,df=df)*tau,col=1)

 1. Normal Linear Model


以下、下記 Lncaster より転載（未テスト）。

model{
for (i in 1:n){
y[i] ~dnorm(mu[i], tau)
mu[i] <- b0 + b1*x[i] + b2*z[i]
}
b0 ~dnorm(0, 1.0E-6)
b1 ~dnorm(0, 1.0E-6)
b2 ~dnorm(0, 1.0E-6)
tau ~dgamma(1.0E-3, 1.0E-3)
}


 2. Heteroscedastic Regression


model{
for (i in 1:n){ y[i] ~dt(mu[i], tau, d)
mu[i] <- b0 + b1*x[i] + b2*z[i]
}
d <- 4
b0 ~dnorm(0, 1.0E-6)
b1 ~dnorm(0, 1.0E-6)
b2 ~dnorm(0, 1.0E-6)
tau ~dgamma(1.0E-3, 1.0E-3)
}

 3. Regression with Autocorrelated Errors


model{
for (i in 2:n){
y[i] ~dnorm(mu[i], tau)
mu[i] <- rho*y[i-1] + b0*(1-rho) + b1*x[i] - rho*b1*x[i-1]
}
b0 ~dnorm(0, 1.0E-6)
b1 ~dnorm(0, 1.0E-6)
tau ~dgamma(1.0E-3, 1.0E-3)
rho ~ dunif(-0.999, 0.999)
}

 4. CES Production Function


model{
for (i in 1:n){
yl[i] ~dnorm(mu[i], tau)
mu[i] <- a - (1/b) * log(d+(1-d)*exp(-b*kl[i]))
}
b <- (1/sig) - 1
sig ~dgamma(2,2)
a ~dnorm(0, 0.001)
d ~dbeta(2,2)
tau ~dgamma(1.0E-3, 1.0E-3)
}

 5. Probit Model


model{
for (i in 1:n){
y[i] ~dbin(p[i], 1)
p[i] <- phi(b0 + b1*x1[i] + b2*x2[i]
}
b0 ~dnorm(0, 1.0E-3)
b1 ~dnorm(0, 1.0E-3)
b2 ~dnorm(0, 1.0E-3)
}

 6 Tobit Model


model{
for (i in 1:n){
ones[i] <- 1
ones[i] ~dbern( p[i] )
term1[i] <- 0.39894*pow(tau, 1/2)*exp(-0.5*tau*pow(y[i] - mu[i], 2))
term2[i] <- phi(-mu[i]*pow(tau, 1/2))
p[i] <- pow(term1[i],ind[i])*pow(term2[i],1-ind[i])/10000
mu[i] <- b0 + b1*x[i]
}
b0 ~dnorm(0, 1.0E-3)
b1 ~dnorm(0, 1.0E-3)
tau ~dgamma(1.0E-3, 1.0E-3)
}

 7. Truncated Normal


model{
for (i in 1:n){
ones[i] <- 1
ones[i] ~dbern( p[i] )
term1[i] <- 0.39894*pow(tau, 1/2)*exp(-0.5*tau*pow(y[i] - mu[i], 2))
term2[i] <- phi(-mu[i]*pow(tau, 1/2))
r[i] <- term1[i]/term2[i]
p[i] <- pow(r[i]-ind[i])/10000
mu[i] <- b0 + b1*x[i]
}
b0 ~dnorm(0, 1.0E-3)
b1 ~dnorm(0, 1.0E-3)
tau ~dgamma(1.0E-3, 1.0E-3)
}

 8. Ordered Probit


model{
for (i in 1:n){
y[i] ~dcat( p[i,] )
p[i, 1] <- phi(r-mu[i])
p[i, 4] <- 1 - phi(r+d+e-mu[i])
p[i, 2] <- phi(r+d-mu[i]) - phi(r-mu[i])
p[i, 3] <- phi(r+d+e-mu[i]) - phi(r+d-mu[i])
mu[i] <- b*x[i]
}
r ~dnorm(0, 1.0E-3)
d ~dgamma(0.1, 0.1)
e ~dgamma(0.1, 0.1)
b ~dnorm(0, 1.0E-3)}
}

 9. Poisson Regression


model{
for(i in 1:n){
y[i] ~ dpois(mu[i])
mu[i] ~ exp(b0 + b1 * x[i])
}
b0 ~ dnorm(0, 0.001)
b1 ~ dnorm(0, 0.001)
}

 10. Heterogeneous Poisson Regression


model{
for(i in 1:n){
y[i] ~ dpois(mu[i])
mu[i] ~ exp(b0 + b1 * x[i]) * v[i]
v[i] ~ dgamma(a, a) 
}
a <- 1/ra
ra ~ dgamma(1, 1)
b0 ~ dnorm(0, 0.001)
b1 ~ dnorm(0, 0.001)
}

 11. Right Censored Weibull Data


model{
for(i in 1:n){
y[i] ~ dweib(alpha, r[i])I(cenind[i], )
mu[i] ~ exp(b0 + b1 * x[i]) 
}
b0 ~ dnorm(0, 0.001)
b1 ~ dnorm(0, 0.001)
alpha ~ dgamma(1, 1)
}

 12. A Censored Heterogeneous Weibull Model


model{
for(i in 1:n){
y[i] ~ dweib(alpha, r[i])I(cenind[i], )
mu[i] ~ v[i] * exp(b0 + b1 * x[i]) 
v[i] ~ dgamma(del, del)
}
del ~ dgamma(1, 1)
alpha ~ dgamma(1, 1)
b0 ~ dnorm(0, 0.001)
b1 ~ dnorm(0, 0.001)
vari <- 1/del
}

 13. An Over-Identified Recursive Equations Model


model{
for (i in 1:n){
y[i,1:2] ~dmnorm(mu[i,], R[,])
mu[i,1] <- b0 + b1*c0*z[i,1] + b1*c2*z[i,2]
mu[i,2] <- c0 + c1*z[i,1] + c2*z[i,2]
}
R[1:2,1:2] ~dwish(Omega[,],2)
b0 ~dnorm(0, 1.0E-4)
b1 ~dnorm(0, 1.0E-4)
c0 ~dnorm(0, 1.0E-4)
c1 ~dnorm(0, 1.0E-4)
c2 ~dnorm(0, 1.0E-4)
}

 14. An Just Identified Simulteneous Equations Model


model{
for (i in 1:n){
y[i,1:2] ~dmnorm(mu[i,], R[,])
mu[i,1] <- ap + bp[1]*zd[i] + bp[2]*zs[i]
mu[i,2] <- aq + bq[1]*zd[i] + bq[2]*zs[i]
}
ap ~dnorm(0, 1.0E-3)
aq ~dnorm(0, 1.0E-3)
for (i in 1:2){
bp[i] ~dnorm(0, 1.0E-3)
bq[i] ~dnorm(0, 1.0E-3)
}
R[1:2,1:2] ~dwish(Omega[,],2)
}

 15. A Panel Data Linear Model


model{
for (i in 1:N){for (t in 1:T){
y[i,t] ~dnorm(mu[i,t],tau)
mu[i,t] <- beta*x[i,t] + alpha[i]
}}
for (i in 1:N){ alpha[i] ~dnorm(alphabar,phi) }
beta  ~dnorm(0, 1.0E-4)
tau ~dgamma(0.1, 0.1)
phi ~dgamma(0.1, 0.1)
}

 16. A Second Order Autoregression


model{
for (i in 3:n){
y[i] ~dnorm(mu[i],tau)
mu[i] <- b0 + b1*y[i-1] + b2*y[i-2]
}
y[1] ~dnorm(0, 1.0E-4)
y[2] ~dnorm(0, 1.0E-4)
b0 ~dnorm(0, 1.0E-4)
b1 ~dnorm(0, 1.0E-4)
b2 ~dnorm(0, 1.0E-4)
tau ~dgamma(1.0E-3, 1.0E-3)
}

 17. Stochastic Volatility


model{
for (i in 1:n){
y[i] ~dnorm(0,p[i])
p[i] <- exp(-h[i])
}
h[1] ~dnorm(muh,qh)
muh <- alpha/(1-rho)
qh <- qv*(1-rho*rho)
for(j in 2:n){
h[j] ~dnorm(mu2[j],pv)
mu2[j] <- alpha + rho*h[j-1]}
alpha ~dnorm(0, 1.0E-4)
tau ~dgamma(1.0E-3, 1.0E-3)
rho ~ dunif(-0.999, 0.999)
}

BUGSコード(2)


 マルコフ切替モデル(単位根検定)


model{
  for(i in 2:n){
    s[i] ~ dcat(pi[s[i-1],])
    mu[i] <- mu0[i,s[i]]
    x[i] ~ dnorm(mu[i],tau[s[i]])
    mu0[i,1] <- phi[2] + (phi[4]+1)*x[i-1]
    mu0[i,2] <- phi[1] + (phi[3]+1)*x[i-1]
  }
  # priors
  phi[1] ~ dnorm(0, 0.001)
  phi[2] <- 0
  phi[3] <- 0
  phi[4] ~ dnorm(0, 0.001)
  for(j in 1:2){
    pi[j,1:2] ~ ddirch(alpha[])
  }
  tau0 ~ dgamma(0.01,0.01)
  omega ~ dgamma(0.01,0.01)
  tau[1] <- tau0
  tau[2] <- tau0*omega
  s[1] ~ dcat(pi0[])
  alpha[1] <- 1
  alpha[2] <- 1
}

関連文献およびリンク



  • Koop, Gary. Bayesian Econometrics. Wiley, 2003.
  • Lncaster, Tony. An Intoroduction to Modern Bayesian Econometrics. Blackwell, 2004.

2冊ともそれほど厚くはありませんが、入門から応用まで広い内容を扱っています。
また、それぞれ扱っている内容や使うプログラミング・ソフトが異なるので、重宝します。使うソフトはKoopがMATLAB、LuncasterはRとWinBUGSです。ただし、実例自体は多くないので別個論文等で補う必要はあるでしょう。

  • 中妻照雄『ファイナンスのためのMCMC法によるベイズ分析』三菱経済研究所, 2003年

  • 和合肇編著『ベイズ計量経済分析』東洋経済新報社, 2005.

基礎から応用まで。基礎編を入門としてつかうには厳しいですが、応用編は論文集になっていて、まだベイズ分析が少数派の現状にあっては、とても参考になります。



Index

講義のページ

著書

経済学

数学・情報

リンク

サイト内検索

キーワード

更新履歴

2022/4/7

2022/3/18

2022/2/21

Copyright (C) 2002 Ryo Hasumi
Powered by FreeStyleWiki3.6.5 with Perl5.008009