3限目

# 時系列回帰モデルを作成
T <- 500
et <- rnorm(T,0,1)
ro <- 0.5
xt <- matrix(0,nrow=1,ncol=T)
xt[1,1] <- et[1]
for(i in 2:T){xt[,i]<-xt[,i-1]*ro+et[i]}
xt.mean <- matrix(0,nrow=1,ncol=T)
xt.mean[1]<-xt[1]
for(i in 2:T){xt.mean[,i]<-sum(xt[,1:i])/i}
par(mfrow=c(3,1))
plot(1:T, xt, type="l", xlab="T", ylab="Xt")
plot(xt[1:T-1], xt[2:T], xlab="Xt", ylab="Xt+1")
plot(1:T, xt.mean,type="l")
# 標本平均の収束過程
T <- 500
et <- matrix(0,nrow=5,ncol=T)
for(i in 1:5){et[1,]<-rnorm(T,0,1)}
ro <- 0.5
xt <- matrix(0,nrow=5,ncol=T)
xt[,1] <- runif(5,-10,10)
for(j in 1:5){
for(i in 2:T){xt[j,i]<-xt[j,i-1]*ro+et[i]}}
xt.mean <- matrix(0,nrow=5,ncol=T)
xt.mean[1]<-xt[1]
for(j in 1:5){
for(i in 2:T){xt.mean[j,i]<-sum(xt[j,1:i])/i}}
plot(xt.mean[1,],type="l", ylim=c(min(xt.mean),max(xt.mean)), xlim=c(1,T), lwd=2)
lines(xt.mean[2,],col="blue", lwd=2)
lines(xt.mean[3,],col="green", lwd=2)
lines(xt.mean[4,],col="red", lwd=2)
lines(xt.mean[5,],col="cyan", lwd=2)
# ギブズ・サンプラー
# 初期値(10, 10)の場合
T<-500; ro<-0.7; y0<-10; mu1<-mu2<-0; sig1<-sig2<-1
y1 <- rep(y0,T); y2<-rep(y0,T)
for(i in 2:T){
y2[i] <- rnorm(1, (mu2+ro*sig2/sig1*(y1[i-1]-mu1)), sqrt(sig2*(1-ro^2)))
y1[i] <- rnorm(1, (mu1+ro*sig1/sig2*(y2[i-1]-mu1)), sqrt(sig1*(1-ro^2)))}
# 初期値(-10, 10)の場合
T<-500; ro<-0.7; y00<--10; mu1<-mu2<-0; sig1<-sig2<-1
y3 <- rep(y0,T); y4<-rep(y00,T)
for(i in 2:T){
y4[i] <- rnorm(1, (mu2+ro*sig2/sig1*(y3[i-1]-mu1)), sqrt(sig2*(1-ro^2)))
y3[i] <- rnorm(1, (mu1+ro*sig1/sig2*(y4[i-1]-mu1)), sqrt(sig1*(1-ro^2)))}
# 二変量正規分布ギブズ・サンプラー
par(mfrow=c(1,2))
plot(y1,y2,type="o",xlim=c(min(y1)-2,max(y1)+2),ylim=c(min(y2)-2,max(y2)+2))
plot(y3,y4,type="o",xlim=c(min(y1)-2,max(y1)+2),ylim=c(min(y2)-2,max(y2)+2))
# 初期値(10, 10)の場合の信頼区間
HPDUPy1<-mean(y1)+qt(0.975,T-2)*sd(y1)
HPDLWy1<-mean(y1)+qt(0.025,T-2)*sd(y1)
HPDUPy1
HPDLWy1
HPDUPy2<-mean(y2)+qt(0.975,T-2)*sd(y2)
HPDLWy2<-mean(y2)+qt(0.025,T-2)*sd(y2)
HPDUPy2
HPDLWy2
# 線形回帰モデルへのギブズ・サンプラーの適用
library(MCMCpack)
data(swiss)
n<-nrow(swiss)
k<-ncol(swiss)
# 通常最小二乗法の適用
sample.lm <- lm(Fertility~.,data=swiss)
summary(sample.lm)
# Gibbs sampling
# burnin=稼働検査期間,mcmc=MCMCの生成回数,b0とB0は事前分布の平均と分散
# c0とd0はτに関するガンマ関数の平均と分散に関するパラメータ
sample.gibbs.post1 <- 
MCMCregress(Fertility~.,data=swiss,
burnin=10000, mcmc=100000,
marginal.likelihood="Chib95", 
b0=0, B0=0.001,c0=0.001, d0=0.001)
summary(sample.gibbs.post1, digit=3)
plot(sample.gibbs.post1)
# Gelman-Rubin統計量
sample.gibbs.post2 <- 
MCMCregress(Fertility~.,data=swiss,burnin=10000, mcmc=100000,
marginal.likelihood="Chib95", b0=1, B0=0.05,c0=0.1, d0=0.1)
sample.gibbs.post.all<-mcmc.list(sample.gibbs.post1,sample.gibbs.post2)
gelman.diag(sample.gibbs.post.all)
# Gewekeの方法
geweke.diag(sample.gibbs.post1)
geweke.diag(sample.gibbs.post2)
# Raftery-Lewisの方法
raftery.diag(sample.gibbs.post1)
# boaパッケージを使う
#library(boa)
# codaパッケージを使う
library(coda)
codamenu()
# ベイズファクターの計算
# model1:最初に推定したモデル
sample.gibbs.post1 <- 
MCMCregress(Fertility~.,data=swiss,
marginal.likelihood="Chib95",
burnin=10000, mcmc=100000, b0=0, B0=0.001,c0=0.001, d0=0.001)
summary(sample.gibbs.post1, digit=3)
# model3:AgricultureとExaminationを除いて推定
sample.gibbs.post3<-
MCMCregress(Fertility~Education+Catholic+Infant.Mortality,data=swiss,
##MCMCregress(Fertility~Education+Catholic+Examination,data=swiss,
marginal.likelihood="Chib95",
burnin=10000, mcmc=100000,
b0=0, B0=0.0015,c0=0.001, d0=0.001)
# ベイズファクター
BayesFactor(sample.gibbs.post1,sample.gibbs.post3)

トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2009-09-14 (月) 01:11:19 (3714d)