贴代码】R语言用metropolis-hastings算法产生4个自由度的t分布的均值,采用如下建议分布:1,N(0,1);2,t(2)

Question

如题：用R语言利用metropolis-hastings算法产生4个自由度的t分布的均值，采用如下建议分布：1、N（0,1）；2、t（2）
最好还能比较一下这两个建议分布
要求贴代码！！！贴代码！！！贴代码！！！

Answer 1

没太懂你说的建议分布是什么意思。比如N(0,1)，是说X^{n+1}=X^{n}+X（X~N(0,1)）还是直接X^{n+1}~N(0,1)？我当是后一种了哈。

#1. N(0,1)
set.seed(2015) #set random seed
N<-1e4 #sample size
x1<-rep(NA,N)
x1[1]<-rnorm(1) #or another initial guess
for(n in 2:N){
    x1[n]<-rnorm(1)
    log.rho<-dt(x1[n],df=4,log=TRUE)-dt(x1[n-1],df=4,log=TRUE)+dnorm(x1[n-1],log=TRUE)-dnorm(x1[n],log=TRUE)
    if(runif(1)>exp(log.rho)) x1[n]<-x1[n-1]
}
mean(x1)

#2. t(2)
N<-1e4
x2<-rep(NA,N)
x2[1]<-rt(1,df=2) #or another initial guess
for(n in 2:N){
    x2[n]<-rt(1,df=2)
    log.rho<-dt(x2[n],df=4,log=TRUE)-dt(x2[n-1],df=4,log=TRUE)+dt(x2[n-1],df=2,log=TRUE)-dt(x2[n],df=2,log=TRUE)
    if(runif(1)>exp(log.rho)) x2[n]<-x2[n-1]
}
mean(x2)

比较的话你就画histogram之类的吧，比如：

hist(x1,freq=FALSE,breaks=40);curve(dt(x,4),-4,4,add=TRUE)
hist(x2,freq=FALSE,breaks=40);curve(dt(x,4),-4,4,add=TRUE)

感觉t(2)要好一些，毕竟proposal不像N(0,1)那样保守。