想在R中实现MCMC的MH算法，找到了MHadaptive包，编写出错。
请问各位大神，在R中是否有其他软件包可以实现MH算法？
求助！
报错为：
Error in li_func(pars, ...) :
unused argument (data = c(13.7168018290208, 2.43176602816527, 10.6765122358815, 15.80358874191, 10.655111449607, -0.195113438451521, 9.6881120593819, 7.72980937745959, 7.42310548938604, 19.7701556113501, 10.7819843821222, 1.83165668335545, 7.60799811728853, 13.8245824558228, 11.0861739773354, 9.87726274886303, 17.8636649924113, 2.8819392145203, 16.949557719703, 14.4423595833012, 23.7879576928322, 14.6814611887133, -7.00768217378476, 1.43957354763575, 0.433094981924043, 14.041879238363, 3.29644328179122,
-5.86551090632578, 0.575169549454536, 6.57002104215313, -7.7726064851087, 11.9942852811918, 5.32637319821304, 5.16529291598935, 17.6351408031294, 16.3909043574436, 21.54023028995, 21.5890958946679, -1.90624858968057, -2.01115126785469, 8.46778703917611, 12.0472208530146, 1.49937266390884, 0.217428318422499, 13.4581775358992, 4.60304879819332, 5.6767169446721, 10.2855181990085, 10.773416663384, 19.5509395698454, 7.7176598401012, -2.61488375161332, 5.70273803412963, 9.676768488
程序如下：
library(MHadaptive)
##各侧面水平数
pnum<-100
inum<-20
##方差分量参数值
var<-c(4,16,64)
##产生多少批数据
iter=1
##储存方差分量、标准误、置信区间和覆盖率
inmc.vc<-array(0,dim=c(iter,length(var)))
inmc.se<-array(0,dim=c(iter,length(var)))
inmc.ci50<-array(0,dim=c(iter,length(var)))
inmc.ci950<-array(0,dim=c(iter,length(var)))
inmc.cover<-array(0,dim=c(iter,length(var)))
#定义似然函数
likelihood <- function(param){
sigma_p<-param[1]
sigma_i<-param[2]
sigma_pi<-param[3]
mu<-array(data = 0, dim = c(p_num,i_num))
log_likelihood<-array(data = 0, dim = c(p_num,i_num))
P<-rnorm(n = p_num, mean = 0, sd = sigma_p)
I<-rnorm(n = i_num, mean = 0, sd = sigma_i)
for (p in 1:p_num) {
for (i in 1:i_num) {
mu[p,i]<-P[p]+I[i]
log_likelihood[p,i]<-dnorm(x = x[p,i], mean = mu[p,i], sd = sigma_pi, log = T)
}
}
sigma_p2<-sigma_p^2
sigma_i2<-sigma_i^2
sigma_pi2<-sigma_pi^2
log_likelihood<-sum(log_likelihood)
return(log_likelihood + prior)
}
#定义先验分布
prior <- function(param){
sigma_p<-param[1]
sigma_i<-param[2]
sigma_pi<-param[3]
prior_sigma_p<-dgamma(sigma_p,2,1/4,log=TRUE)
prior_sigma_i<-dgamma(sigma_i,2,1/16,log=TRUE)
prior_sigma_pi<-dgamma(sigma_pi,2,1/64,log=TRUE)
return(prior_sigma_p + prior_sigma_i + prior_sigma_pi)
}
##模拟正态分布的数据
file_path1<-"C:/Users/Administrator/Desktop/infmc/pi/data/data_"
for (k in 1:iter){
pz<-rnorm(pnum,0,(var[1])^0.5)
iz<-rnorm(inum,0,(var[2])^0.5)
piz<-array(rnorm((pnum*inum),0,(var[3])^0.5),dim=c(pnum,inum))
d<-array(0,dim=c(pnum,inum))
for(p in 1:pnum){
for(i in 1:inum){
d[p,i]<-0+pz[p]+iz[i]+piz[p,i]
}
}
write.csv(d,file=paste(file_path1,k,".csv",sep = ""),row.names = F )
##MCMC方法（有信息先验）
mcmc_r<-Metro_Hastings(li_func=likelihood,pars=c(0.01,0.01,0.01),
par_names=c('varp','vari','vare'),iterations = 10000, burn_in = 5000,
adapt_par = c(5000, 20, 0.5, 0.75),data=d)
result<-BCI(mcmc_r, interval = c(0.05, 0.95))
plotMH(mcmc_r)
inmc.vc[k,]<-result[1:3,3]
inmc.ci50[k,]<-result[1:3,1]
inmc.ci950[k,]<-result[1:3,2]
inmc.cover[k,]<-ifelse(inmc.ci50[k,]<var&var<inmc.ci950[k,],1,0)
}
