Go to Qiita Advent Calendar 2024 Top
LoginSignup
2
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 3 years have passed since last update.

@kozakai-ryouta(亮太 小酒井)

ベイズ計量経済分析（マルコフ連鎖モンテカルロ法とその応用） 和合肇先生

Last updated at Posted at 2021-08-27

#正規分布（ギブスサンプラー） p.45

library(MCMCpack)

n0=10;n=50;n1=n0+n;S0=1

x=rnorm(n,25,6)

mu0=1;sig0=10

ite=10^4

for(l in 1:ite){
  
if(l==1){
  
sig1=sqrt(1/(1/(sig0^2)+n/(sig0^2)))

mu1=(sig1^2)*(mu0/(sig0^2)+n*mean(x)/(sig0^2))

mu=mu1;sig=sig1  
  
}else{
  
sig1=sqrt(1/(1/(sig0^2)+n/(sig^2)))

mu1=(sig1^2)*(mu0/(sig0^2)+n*mean(x)/(sig^2))

mu=mu1;sig=sig1  
  
}  

mu=rnorm(1,mu,sig^2)

sig <- sqrt(rinvgamma(1,shape=n1/2,scale=(n0*S0+sum((x-mu)^2))/2))

print(paste0(mu,":",sig))

}



#構造変化のある回帰model　p.49

library(MASS);library(MCMCpack)

data(iris)

iris=iris[1:100,]

X=as.matrix(iris[,!(colnames(iris) %in% "Species")])/100

Y=c(rep(10,50),rep(15,50))/1

k=5

n=nrow(X);n0=2;m0=2;S0=1;V0=1

b0=rep(1,ncol(X));a0=rep(1,ncol(X))

beta=rep(1,ncol(X))/1000;alpha=rep(1,ncol(X))/1000

B0=diag(1,ncol(X));A0=diag(1,ncol(X))

pthi1=1;pthi2=1




ite=100

k_vec=c()

for(l in 1:ite){

X1=X[1:k,];X2=X[-c(1:k),]

mat1=array(0,dim=c(ncol(X),ncol(X)))

for(i in 1:nrow(X1)){

mat1=mat1+t(t(c(X1[i,])))%*%t(c(X1[i,]))  

}

mat2=array(0,dim=c(ncol(X),ncol(X)))

for(i in 1:nrow(X2)){

mat2=mat2+t(t(c(X2[i,])))%*%t(c(X2[i,]))  

}

A1=solve(solve(A0)+pthi1*mat1)

B1=solve(solve(B0)+pthi2*mat2)

a1=A1%*%(solve(A0)%*%a0+pthi1*t(X1)%*%Y[1:k])

b1=B1%*%(solve(B0)%*%b0+pthi2*t(X2)%*%Y[-c(1:k)])

beta=(mvrnorm(1,b1,B1))

alpha=(mvrnorm(1,a1,A1))

n1=n0+k;m1=m0+(n-k)

pthi1=mean(rgamma(1,n1/2,(n0*S0+sum((Y[1:k]-X1%*%alpha)^2))/2))

pthi2=mean(rgamma(1,m1/2,(m0*V0+sum((Y[-c(1:k)]-X2%*%beta)^2))/2))

w=c()

for(j in 1:n){

val=(pthi1^(j/2))*(pthi2^((n-j)/2))*exp(-(pthi1*sum((Y[1:j]-X[1:j,]%*%alpha)^2)+pthi2*sum((Y[-c(1:j)]-X[-c(1:j),]%*%beta)^2))/2)  

w=c(w,val)

}

p=w/sum(w);

#k=c(1:n)[p==max(p)]

k=sample(1:n,1,prob=p)

k=ifelse(k<5,5,ifelse(k>n-5,n-5,k))

k_vec=c(k_vec,k)

cost=-(pthi1*sum((Y[1:k]-X[1:k,]%*%alpha)^2)+pthi2*sum((Y[-c(1:k)]-X[-c(1:k),]%*%beta)^2))/2

print(cost)


}

plot(k_vec)



#t分布回帰model　p.56

library(MASS);library(MCMCpack)

data=data.frame(num=1:20,y=c(0.20, 0.10, 0.49, 0.26, 0.92, 0.95, 0.18, 0.19, 0.44, 0.79, 0.61, 0.41, 0.49, 0.34, 0.62, 0.99, 0.38, 0.22,0.71,0.40),x1=c(84,87,86,85,82,87,92,94,88,84.9,78,90,88,87,82,84,86,83,85.2,82),x2=c(61,55.5,57,57,50,50,66.5,65,60.5,49.5,49.5,61,59.5,58.4,53.5,54,60,58.8,54,56),x3=c(24.5,24,23,22.5,22,24,25,26,25,24,23,21,26,25.5,24,23,24,24,24,22),sign=0)

X=as.matrix(cbind(rep(1,nrow(data)),data[,colnames(data) %in% c("x1","x2")]))/10

Y=data$y*1

n=nrow(X);n0=3

b0=rep(1,ncol(X))

beta=rep(1,ncol(X))

B0=diag(1,ncol(X))

lam=diag(1,nrow(X))

pthi=1

v0=1;n1=n0+n

ite=1000

for(l in 1:ite){

B1=solve(solve(B0)+pthi*t(X)%*%lam%*%X)

b1=B1%*%(solve(B0)%*%b0+pthi*t(X)%*%lam%*%Y)

beta=mvrnorm(1,b1,B1)

pthi=rgamma(1,shape=n1/2,scale = (n0*S0+t(Y-X%*%beta)%*%lam%*%(Y-X%*%beta))/2)

lam_vec=c()

for(i in 1:nrow(X)){

lam_vec=c(lam_vec,rgamma(1,shape=(v0+1)/2,scale=(v0+pthi*(Y-X%*%beta)[i]^2)/2))

}

lam=diag(lam_vec)

print(sum((Y-X%*%beta)^2))

}



#yiがベクトルの回帰model　p.58

library(MASS);library(MCMCpack)

data(anscombe)

X=as.matrix(anscombe[,c(6:8)])

Y=as.matrix(anscombe[,c(3:5)])

n=nrow(X);n0=3

b0=matrix(rep(1,ncol(Y)*ncol(X)),ncol=ncol(Y))

beta=matrix(rep(1,ncol(Y)*ncol(X)),ncol=ncol(Y))

B0=diag(1,ncol(X))

pthi=diag(1,nrow(X))

n1=n0+n

R0=diag(1,nrow(X))

ite=10^5

for(l in 1:ite){
  
beta_pre=beta;pthi_pre=pthi

B1=solve(solve(B0)+t(X)%*%pthi%*%X)

b1=B1%*%(solve(B0)%*%b0+t(X)%*%pthi%*%Y)

for(i in 1:ncol(beta)){

beta[,i]=mvrnorm(1,b1[,i],B1)

}

mat=array(0,dim=dim(R0))

for(i in 1:ncol(beta)){
  
mat=mat+t(t((Y-X%*%beta)[,i]))%*%(t((Y-X%*%beta)[,i]))
  
}

pthi=rwish(n1,solve(solve(R0)+mat))

print(sum(abs(beta_pre-beta))+sum(abs(pthi_pre-pthi)))

}



#MH-AR(1) p.67

library(MCMCpack)

data("AirPassengers")

Y=c(AirPassengers)

Y0=Y[1]*0.8

n0=3;n=length(Y);n1=n0+n+1

S0=1

pthi=0.0001

sig=1

ite=1000

for(l in 1:ite){
  
sig_pre=sig

sig=sqrt(rinvgamma(1,shape=n1/2,scale=(n0*S0+(1-pthi^2)*Y0^2+sum((Y-c(Y0,Y)[1:n]*pthi)^2))/2))

pthi_hat=sum(Y*c(Y0,Y)[1:n])/sum(Y[1:(n-1)]^2)

mu=pthi_hat;sigma2=(sig^2)/sum(Y[1:(n-1)]^2)

aa=1;bb=-1;pthi_pre=pthi

pthi=qnorm(runif(1,pnorm(bb,mu,sigma2,1),pnorm(aa,mu,sigma2)),mu,sigma2)

a=min(c(1,sqrt(1-pthi^2)/sqrt(1-pthi_pre^2)))

u=runif(1,0,1)

pthi=ifelse(u<a,pthi,pthi_pre)

print(sum(abs(pthi-pthi_pre))+abs(sig-sig_pre))

}



#ロジット・モデル(MHアルゴリズム) p.69

library(MASS)

data=data.frame(num=1:20,y=c(0.20, 0.10, 0.49, 0.26, 0.92, 0.95, 0.18, 0.19, 0.44, 0.79, 0.61, 0.41, 0.49, 0.34, 0.62, 0.99, 0.38, 0.22,0.71,0.40),x1=c(84,87,86,85,82,87,92,94,88,84.9,78,90,88,87,82,84,86,83,85.2,82),x2=c(61,55.5,57,57,50,50,66.5,65,60.5,49.5,49.5,61,59.5,58.4,53.5,54,60,58.8,54,56),x3=c(24.5,24,23,22.5,22,24,25,26,25,24,23,21,26,25.5,24,23,24,24,24,22),sign=0)

for(j in 1:nrow(data)){

data$sign[j]=ifelse(data$y[j]>0.49,1,0)

}

X=as.matrix(cbind(rep(1,nrow(data)),data[,colnames(data) %in% c("x1","x2")]/10))

Y=data$sign

beta=rep(1,ncol(X))*0.001

p=1-1/(1+exp(-X%*%beta))

b0=t(X)%*%log(p/(1-p))

B0=diag(1,ncol(X))

eta=0.001

ite=1000

for(l in 1:ite){

p=1-1/(1+exp(-X%*%beta))

beta_hat=mvrnorm(1,beta,B0)

beta_pre=beta

#beta更新用

g=sum(Y*log(ifelse(p>0,p,1))+(1-Y)*log(ifelse(1-p>0,1-p,1)))-t(c(beta-b0))%*%solve(B0)%*%t(t(c(beta-b0)))/2

dg=t(X)%*%(Y-p)-(solve(B0)+t(solve(B0)))%*%(beta-b0)/2

ddg=-t(X)%*%diag(c(p))%*%(diag(rep(1,length(Y)))-diag(c(p)))%*%X-(solve(B0)+t(solve(B0)))/2

h=g+sum(dg*(beta-beta_hat))+t(c(beta-beta_hat))%*%ddg%*%t(t(c(beta-beta_hat)))/2

m=beta_hat-eta*solve(ddg)%*%dg;g1=g

beta=mvrnorm(1,m,-solve(ddg))

p=1-1/(1+exp(-X%*%beta))

g2=sum(Y*log(ifelse(p>0,p,1))+(1-Y)*log(ifelse(1-p>0,1-p,1)))-t(c(beta-b0))%*%solve(B0)%*%t(t(c(beta-b0)))/2

#ｈ(beta)計算用

p=1-1/(1+exp(-X%*%beta_pre))

ddg1=-t(X)%*%diag(c(p))%*%(diag(rep(1,length(Y)))-diag(c(p)))%*%X-(solve(B0)+t(solve(B0)))/2

dg1=t(X)%*%(Y-p)-(solve(B0)+t(solve(B0)))%*%(beta_pre-b0)/2

p=1-1/(1+exp(-X%*%beta))

ddg2=-t(X)%*%diag(c(p))%*%(diag(rep(1,length(Y)))-diag(c(p)))%*%X-(solve(B0)+t(solve(B0)))/2

dg2=t(X)%*%(Y-p)-(solve(B0)+t(solve(B0)))%*%(beta-b0)/2

h1=g1+sum(dg1*(beta_pre-beta_hat))+t(c(beta_pre-beta_hat))%*%ddg1%*%t(t(c(beta_pre-beta_hat)))/2

h2=g2+sum(dg2*(beta-beta_hat))+t(c(beta-beta_hat))%*%ddg2%*%t(t(c(beta-beta_hat)))/2

u=runif(1,0,1)

alpha=min(c(exp(g2-g1-(h2-h1)),1))

if(u<alpha){beta=beta}else(beta=beta_pre)

print(beta)

}

cbind(Y,p)




rm(list=ls())

#ロジット・モデル(MHアルゴリズム) p.69

library(MASS)

data=data.frame(num=1:20,y=c(0.20, 0.10, 0.49, 0.26, 0.92, 0.95, 0.18, 0.19, 0.44, 0.79, 0.61, 0.41, 0.49, 0.34, 0.62, 0.99, 0.38, 0.22,0.71,0.40),x1=c(84,87,86,85,82,87,92,94,88,84.9,78,90,88,87,82,84,86,83,85.2,82),x2=c(61,55.5,57,57,50,50,66.5,65,60.5,49.5,49.5,61,59.5,58.4,53.5,54,60,58.8,54,56),x3=c(24.5,24,23,22.5,22,24,25,26,25,24,23,21,26,25.5,24,23,24,24,24,22),sign=0)

for(j in 1:nrow(data)){

data$sign[j]=ifelse(data$y[j]>0.49,1,0)

}

X=as.matrix(cbind(rep(1,nrow(data)),data[,colnames(data) %in% c("x1","x2")]/10))

Y=data$sign

beta=rep(1,ncol(X))*0.001

p=1-1/(1+exp(-X%*%beta))

b0=t(X)%*%log(p/(1-p))

B0=diag(1,ncol(X))

eta=1

ite=1000

pridis=function(b){

prob=1-1/(1+exp(-X%*%b))

loglik=sum(Y*log(ifelse(prob>0,prob,1))+(1-Y)*log(ifelse(1-prob>0,1-prob,1)))

dis=loglik-t(c(b-b0))%*%solve(B0)%*%t(t(c(b-b0)))

return(dis)

}


for(l in 1:ite){

p=1-1/(1+exp(-X%*%beta))

samp=mvrnorm(100,beta,diag(1,length(beta)))

loglik_vec=c()

for(j in 1:100){

loglik_vec=c(loglik_vec,pridis(c(samp[j,])))

}

ps=abs(loglik_vec)/sum(abs(loglik_vec))

beta_hat=samp[sample(c(1:100),1,prob=ps),]

#beta_hat=mvrnorm(1,beta,B0)

beta_pre=beta

#beta更新用

g=sum(Y*log(ifelse(p>0,p,1))+(1-Y)*log(ifelse(1-p>0,1-p,1)))-t(c(beta-b0))%*%solve(B0)%*%t(t(c(beta-b0)))/2

dg=t(X)%*%(Y-p)-(solve(B0)+t(solve(B0)))%*%(beta-b0)/2

ddg=-t(X)%*%diag(c(p))%*%(diag(rep(1,length(Y)))-diag(c(p)))%*%X-(solve(B0)+t(solve(B0)))/2

h=g+sum(dg*(beta-beta_hat))+t(c(beta-beta_hat))%*%ddg%*%t(t(c(beta-beta_hat)))/2

m=beta_hat-eta*solve(ddg)%*%dg;g1=g

beta=mvrnorm(1,m,-solve(ddg))

p=1-1/(1+exp(-X%*%beta))

g2=sum(Y*log(ifelse(p>0,p,1))+(1-Y)*log(ifelse(1-p>0,1-p,1)))-t(c(beta-b0))%*%solve(B0)%*%t(t(c(beta-b0)))/2

#ｈ(beta)計算用

p=1-1/(1+exp(-X%*%beta_pre))

ddg1=-t(X)%*%diag(c(p))%*%(diag(rep(1,length(Y)))-diag(c(p)))%*%X-(solve(B0)+t(solve(B0)))/2

dg1=t(X)%*%(Y-p)-(solve(B0)+t(solve(B0)))%*%(beta_pre-b0)/2

p=1-1/(1+exp(-X%*%beta))

ddg2=-t(X)%*%diag(c(p))%*%(diag(rep(1,length(Y)))-diag(c(p)))%*%X-(solve(B0)+t(solve(B0)))/2

dg2=t(X)%*%(Y-p)-(solve(B0)+t(solve(B0)))%*%(beta-b0)/2

h1=g1+sum(dg1*(beta_pre-beta_hat))+t(c(beta_pre-beta_hat))%*%ddg1%*%t(t(c(beta_pre-beta_hat)))/2

h2=g2+sum(dg2*(beta-beta_hat))+t(c(beta-beta_hat))%*%ddg2%*%t(t(c(beta-beta_hat)))/2

u=runif(1,0,1)

alpha=min(c(exp(g2-g1-(h2-h1)),1))

if(u<alpha){beta=beta}else(beta=beta_pre)

print(sum(abs(Y-ifelse(p>0.5,1,0)))/length(Y))

}

cbind(Y,ifelse(p>0.5,1,0))




#ロジット・モデル(AR-MHアルゴリズム) p.70

library(MASS)

data=data.frame(num=1:20,y=c(0.20, 0.10, 0.49, 0.26, 0.92, 0.95, 0.18, 0.19, 0.44, 0.79, 0.61, 0.41, 0.49, 0.34, 0.62, 0.99, 0.38, 0.22,0.71,0.40),x1=c(84,87,86,85,82,87,92,94,88,84.9,78,90,88,87,82,84,86,83,85.2,82),x2=c(61,55.5,57,57,50,50,66.5,65,60.5,49.5,49.5,61,59.5,58.4,53.5,54,60,58.8,54,56),x3=c(24.5,24,23,22.5,22,24,25,26,25,24,23,21,26,25.5,24,23,24,24,24,22),sign=0)

for(j in 1:nrow(data)){

data$sign[j]=ifelse(data$y[j]>0.49,1,0)

}

X=as.matrix(cbind(rep(1,nrow(data)),data[,colnames(data) %in% c("x1","x2")]/10))

Y=data$sign


b0=rep(1,ncol(X))

beta=rep(1,ncol(X))*0.01

B0=diag(1,ncol(X))*0.001


ite=10

for(l in 1:ite){

p=1-1/(1+exp(-X%*%beta))

beta_hat=mvrnorm(1,b0,B0)

#beta更新用

dg=t(X)%*%(Y-p)+(solve(B0)+t(solve(B0)))%*%(beta-b0)/2

ddg=array(0,dim=c(ncol(X),ncol(X)))

ddg=-t(X)%*%diag(c(p))%*%(diag(rep(1,length(Y)))-diag(c(p)))%*%X+(solve(B0)+t(solve(B0)))/2

m=beta+solve(ddg)%*%dg

g1=sum(Y*(-X%*%beta)-log(1+exp(-X%*%beta)))-t(c(beta-b0))%*%solve(B0)%*%t(t(c(beta-b0)))/2

beta_pre=beta

beta=mvrnorm(1,m,ddg)

g2=sum(Y*(-X%*%beta)-log(1+exp(-X%*%beta)))-t(c(beta-b0))%*%solve(B0)%*%t(t(c(beta-b0)))/2

#ｈ(beta)計算用

p=1-1/(1+exp(-X%*%beta_pre))

ddg1=-t(X)%*%diag(c(p))%*%(diag(rep(1,length(Y)))-diag(c(p)))%*%X+(solve(B0)+t(solve(B0)))/2

dg1=t(X)%*%(Y-p)+(solve(B0)+t(solve(B0)))%*%(beta_pre-b0)/2

p=1-1/(1+exp(-X%*%beta))

ddg2=-t(X)%*%diag(c(p))%*%(diag(rep(1,length(Y)))-diag(c(p)))%*%X+(solve(B0)+t(solve(B0)))/2

dg2=t(X)%*%(Y-p)+(solve(B0)+t(solve(B0)))%*%(beta-b0)/2

h1=g1+sum(dg1*(beta_pre-beta_hat))+t(c(beta_pre-beta_hat))%*%ddg1%*%t(t(c(beta_pre-beta_hat)))/2

h2=g2+sum(dg2*(beta-beta_hat))+t(c(beta-beta_hat))%*%ddg2%*%t(t(c(beta-beta_hat)))/2

u=runif(1,0,1)

alpha=min(min(c(1,exp(h1-g1)))/min(c(1,h2-g2)),1)

beta=ifelse(u<alpha,beta,beta_pre)

print(sum(abs(beta-beta_pre)))

}
2
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
2
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?