独立抽样（MCMC方法）

独立抽样（一）

# simulations
wais<-read.csv("wais.csv",head=T)
y<-wais[,2]
x<-wais[,1]
m<-55000
mu.beta<-c(0,0)
sigma.beta<-c(100,100)
prop.s<-c(.1,.1)
beta<-matrix(nrow=m,ncol=2)
acc.prob<-0
current.beta<-c(0,0)
for(t in 1:m){
  prop.beta<-rnorm(2,current.beta,prop.s)
 cur.eta<-current.beta[1]+current.beta[2]*x
  prop.eta<-prop.beta[1]+prop.beta[2]*x
 loga<-(sum(y*prop.eta-log(1+exp(prop.eta))))-sum(y*cur.eta-log(1+exp(cur.eta)))+sum(dnorm(prop.beta,mu.beta,prop.s,log=T))-sum(dnorm(current.beta,mu.beta,prop.s,log=T))
  u<-runif(1)
  u<-log(u)
  if(u<=loga){
   current.beta<-prop.beta
   acc.prob<-acc.prob+1
  }
  beta[t,]<-current.beta
}
    #convergence diagnostic plot
erg.mean<-function(x){
  n<-length(x)
  result<-cumsum(x)/cumsum(rep(1,n))
}
burnin<-15000
idx<-seq(1,m,50)
idx2<-seq(burnin+1,m)
par(mfrow=c(2,2))
plot(idx,beta[idx,1],type="l",xlab="Iterations",ylab="Values ofbeta0")
plot(idx,beta[idx,2],type="l",xlab="Iterations",ylab="Values ofbeta1")
ergbeta0<-erg.mean(beta[,1])
ergbeta02<-erg.mean(beta[idx2,1])
ylims0<-range(c(ergbeta0,ergbeta02))
ergbeta1<-erg.mean(beta[,2])
ergbeta12<-erg.mean(beta[idx2,2])
ylims1<-range(c(ergbeta1,ergbeta12))
plot(idx,ergbeta0[idx],type="l",ylab='Values ofbeta0',xlab='Iterations',main='(c)Ergodic Mean Plot ofbeta0',ylim=ylims0)
lines(idx2,ergbeta02[idx2-burnin],col=2,lty=2)
plot(idx,ergbeta1[idx],type="l",ylab='Values ofbeta0',xlab='Iterations',main='(d) Ergodic Mean Plot ofbeta0',ylim=ylims1)
lines(idx2,ergbeta12[idx2-burnin],col=2,lty=2)
par(mfrow=c(1,1))
apply(beta[(burnin+1):m,],2,mean)
apply(beta[(burnin+1):m,],2,sd)
cor(beta[(burnin+1):m,1],beta[(burnin+1):m,2])                                           

独立抽样（二）

m<-5000
xt<-numeric(m)
a<-5
b<-2
p<-.2
n<-30
mu<-c(0,5)
sigma<-c(1,1)
i<-sample(1:2,size=n,replace=T,prob=c(p,1-p))
x<-rnorm(n,mu[i],sigma[i])
u<-runif(m)
y<-rbeta(m,a,b)
xt[1]<-.5
for(i in 2:m){
 fy<-y[i]*dnorm(x,mu[1],sigma[1])+(1-y[i])*dnorm(x,mu[2],sigma[2])
 fx<-xt[i-1]*dnorm(x,mu[1],sigma[1])+(1-xt[i-1])*dnorm(x,mu[2],sigma[2])
 r<-prod(fy/fx)*(xt[i-1]^(a-1)*(1-xt[i-1])^(b-1))/(y[i]^(a-1)*(1-y[i])^(b-1))
  if(u[i]<=r) xt[i]<-y[i]
  else xt[i]<-xt[i-1]
}
plot(xt,type="l",ylab="p")
hist(xt[101:m],main="",xlab="p",prob=T)
print(mean(xt[101:m]))