adaptMetropGibbs(spBayes)
adaptMetropGibbs()所属R语言包：spBayes
Adaptive Metropolis within Gibbs algorithm
Gibbs算法的自适应大都市内
译者：生物统计家园网 机器人LoveR
描述----------Description----------
Markov chain Monte Carlo for continuous random vector using an adaptive Metropolis within Gibbs algorithm.
马尔可夫链蒙特卡罗连续随机向量内使用自适应都会Gibbs算法。
用法----------Usage----------
adaptMetropGibbs(ltd, starting, tuning=1, accept.rate=0.44,
batch = 1, batch.length=25, report=100,
verbose=TRUE, ...)
参数----------Arguments----------
参数：ltd
an R function that evaluates the log target density of the desired equilibrium distribution of the Markov chain. First argument is the starting value vector of the Markov chain. Pass variables used in the ltd via the ... argument of aMetropGibbs.
一个R函数计算结果的log目标马尔可夫链所需的均衡分布密度。第一个参数是初始值向量的马氏链。通过使用的变量在ltd通过...参数aMetropGibbs。
参数：starting
a real vector of parameter starting values.
实数向量的参数的初始值。
参数：tuning
a scalar or vector of initial Metropolis tuning values. The vector must be of length(starting). If a scalar is passed then it is expanded to length(starting).
一个标量或矢量的初步大都会的调整值。向量必须是length(starting)。如果一个标量被传递，则它被扩展为length(starting)。
参数：accept.rate
a scalar or vector of target Metropolis acceptance rates. The vector must be of length(starting). If a scalar is passed then it is expanded to length(starting).
一个标量或矢量的目标都会接受率。向量必须是length(starting)。如果一个标量被传递，则它被扩展为length(starting)。
参数：batch
the number of batches.
的批数。
参数：batch.length
the number of sampler iterations in each batch.
每个批处理中的采样迭代的数目。
参数：report
the number of batches between acceptance rate reports.
批次之间的录取率报告。
参数：verbose
if TRUE, progress of the sampler is printed to the screen. Otherwise, nothing is printed to the screen.
如果TRUE，进步的取样打印到屏幕上。否则，一切都被打印到屏幕上。
参数：...
currently no additional arguments.
目前没有任何额外的参数。
值----------Value----------
A list with the following tags:
下面的标签列表：
参数：p.samples
a coda object of posterior samples for the parameters.
coda对象后样品的参数。
参数：acceptance
the Metropolis acceptance rate at the end of each batch.
每一批次结束时都会接受率。
参数：ltd
ltd
ltd
参数：accept.rate
accept.rate
accept.rate
参数：batch
batch
batch
参数：batch.length
batch.length
batch.length
参数：proc.time
the elapsed CPU and wall time (in seconds).
CPU和墙壁的经过的时间（以秒为单位）。
注意----------Note----------
This function is a rework of Rosenthal (2007) with some added
这个函数是一个罗森塔尔（2007年）增加了一些返工
（作者）----------Author(s)----------
Andrew O. Finley <a href="mailto:finleya@msu.edu">finleya@msu.edu</a>, <br>
Sudipto Banerjee <a href="mailto:sudiptob@biostat.umn.edu">sudiptob@biostat.umn.edu</a>
参考文献----------References----------
Statistics and Data Analysis. 51:5467-5470.
实例----------Examples----------
## Not run: [＃不运行：]
###########################[＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃]
##Fit a spatial regression[＃适合的空间回归]
###########################[＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃]
set.seed(1)
n <- 50
x <- runif(n, 0, 100)
y <- runif(n, 0, 100)
D <- as.matrix(dist(cbind(x, y)))
phi <- 3/50
sigmasq <- 50
tausq <- 20
mu <- 150
s <- (sigmasq*exp(-phi*D))
w <- mvrnorm(1, rep(0, n), s)
Y <- mvrnorm(1, rep(mu, n) + w, tausq*diag(n))
X <- as.matrix(rep(1, length(Y)))
##Priors[＃先验的]
##IG sigma^2 and tau^2[IG西格玛＃^ 2和tau ^ 2]
a.sig <- 2
b.sig <- 100
a.tau <- 2
b.tau <- 100
##Unif phi[＃UNIF披]
a.phi <- 3/100
b.phi <- 3/1
##Functions used to transform phi to continuous support.[＃函数用于将披继续支持。]
logit <- function(theta, a, b){log((theta-a)/(b-theta))}
logit.inv <- function(z, a, b){b-(b-a)/(1+exp(z))}
##Metrop. target[＃Metrop。目标]
target <- function(theta){
mu.cand <- theta[1]
sigmasq.cand <- exp(theta[2])
tausq.cand <- exp(theta[3])
phi.cand <- logit.inv(theta[4], a.phi, b.phi)
Sigma <- sigmasq.cand*exp(-phi.cand*D)+tausq.cand*diag(n)
SigmaInv <- chol2inv(chol(Sigma))
logDetSigma <- determinant(Sigma, log=TRUE)$modulus[1]
out <- (
##Priors[＃先验的]
-(a.sig+1)*log(sigmasq.cand) - b.sig/sigmasq.cand
-(a.tau+1)*log(tausq.cand) - b.tau/tausq.cand
##Jacobians[＃雅可比矩阵]
+log(sigmasq.cand) + log(tausq.cand)
+log(phi.cand - a.phi) + log(b.phi -phi.cand)
##Likelihood[＃似然]
-0.5*logDetSigma-0.5*(t(Y-X%*%mu.cand)%*%SigmaInv%*%(Y-X%*%mu.cand))
)
return(out)
}
##Run a couple chains[＃执行一对夫妇链]
n.batch <- 500
batch.length <- 25
inits <- c(0, log(1), log(1), logit(3/10, a.phi, b.phi))
chain.1 <- adaptMetropGibbs(ltd=target, starting=inits,
batch=n.batch, batch.length=batch.length, report=100)
inits <- c(500, log(100), log(100), logit(3/90, a.phi, b.phi))
chain.2 <- adaptMetropGibbs(ltd=target, starting=inits,
batch=n.batch, batch.length=batch.length, report=100)
##Check out acceptance rate just for fun[＃检查验收率只是为了好玩]
plot(mcmc.list(mcmc(chain.1$acceptance), mcmc(chain.2$acceptance)))
##Back transform[＃返回改造]
chain.1$p.samples[,2] <- exp(chain.1$p.samples[,2])
chain.1$p.samples[,3] <- exp(chain.1$p.samples[,3])
chain.1$p.samples[,4] <- 3/logit.inv(chain.1$p.samples[,4], a.phi, b.phi)
chain.2$p.samples[,2] <- exp(chain.2$p.samples[,2])
chain.2$p.samples[,3] <- exp(chain.2$p.samples[,3])
chain.2$p.samples[,4] <- 3/logit.inv(chain.2$p.samples[,4], a.phi, b.phi)
par.names <- c("mu", "sigma.sq", "tau.sq", "effective range (-log(0.05)/phi)")
colnames(chain.1$p.samples) <- par.names
colnames(chain.2$p.samples) <- par.names
##Discard burn.in and plot and do some convergence diagnostics[，丢弃burn.in和图，并做一些收敛诊断]
chains <- mcmc.list(mcmc(chain.1$p.samples), mcmc(chain.2$p.samples))
plot(window(chains, start=as.integer(0.5*n.batch*batch.length)))
gelman.diag(chains)
##########################[＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃]
##Example of fitting a[＃示例的拟合]
##a non-linear model[＃非线性模型]
##########################[＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃]
##Example of fitting a non-linear model[＃拟合非线性模型的实施例]
set.seed(1)
########################################################[＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃ ＃＃＃＃＃]
##Simulate some data.[＃模拟一些数据。]
########################################################[＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃ ＃＃＃＃＃]
a <- 0.1 #-Inf < a < Inf[-INF <A <天道酬勤]
b <- 0.1 #b > 0[B> 0]
c <- 0.2 #c > 0[C> 0，]
tau.sq <- 0.1 #tau.sq > 0[tau.sq> 0]
fn <- function(a,b,c,x){
a+b*exp(x/c)
}
n <- 200
x <- seq(0,1,0.01)
y <- rnorm(length(x), fn(a,b,c,x), sqrt(tau.sq))
##check out your data[＃检查你的数据]
plot(x, y)
########################################################[＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃ ＃＃＃＃＃]
##The log target density[＃log目标密度]
########################################################[＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃ ＃＃＃＃＃]
##Define the log target density used in the Metrop.[＃定义log目标密度在Metrop。]
ltd <- function(theta){
##extract and transform as needed[＃提取和转换，根据需要]
a <- theta[1]
b <- exp(theta[2])
c <- exp(theta[3])
tau.sq <- exp(theta[4])
y.hat <- fn(a, b, c, x)
##likelihood[＃可能性]
logl <- sum(dnorm(y, y.hat, sqrt(tau.sq), log=TRUE))
##priors IG on tau.sq and normal on a and transformed b, c, d[＃先验IG tau.sq和正常的一个转化的B，C，D]
logl <- (logl
-(IG.a+1)*log(tau.sq)-IG.b/tau.sq
+sum(dnorm(theta[1:3], N.mu, N.v, log=TRUE))
)
##Jacobian adjustment for tau.sq[＃雅可比为tau.sq调整]
logl <- logl+log(tau.sq)
return(logl)
}
########################################################[＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃ ＃＃＃＃＃]
##The rest[＃截断]
########################################################[＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃ ＃＃＃＃＃]
##Priors[＃先验的]
IG.a <- 2
IG.b <- 0.01
N.mu <- 0
N.v <- 10
theta.tuning <- c(0.01, 0.01, 0.005, 0.01)
##Run three chains with different starting values[＃执行三个链，用不同的初始值]
n.batch <- 1000
batch.length <- 25
theta.starting <- c(0, log(0.01), log(0.6), log(0.01))
run.1 <- adaptMetropGibbs(ltd=ltd, starting=theta.starting, tuning=theta.tuning,
batch=n.batch, batch.length=batch.length, report=100)
theta.starting <- c(1.5, log(0.05), log(0.5), log(0.05))
run.2 <- adaptMetropGibbs(ltd=ltd, starting=theta.starting, tuning=theta.tuning,
batch=n.batch, batch.length=batch.length, report=100)
theta.starting <- c(-1.5, log(0.1), log(0.4), log(0.1))
run.3 <- adaptMetropGibbs(ltd=ltd, starting=theta.starting, tuning=theta.tuning,
batch=n.batch, batch.length=batch.length, report=100)
##Back transform[＃返回改造]
samples.1 <- cbind(run.1$p.samples[,1], exp(run.1$p.samples[,2:4]))
samples.2 <- cbind(run.2$p.samples[,1], exp(run.2$p.samples[,2:4]))
samples.3 <- cbind(run.3$p.samples[,1], exp(run.3$p.samples[,2:4]))
samples <- mcmc.list(mcmc(samples.1), mcmc(samples.2), mcmc(samples.3))
##Summary [＃摘要]
plot(samples, density=FALSE)
gelman.plot(samples)
burn.in <- 5000
fn.pred <- function(theta,x){
a <- theta[1]
b <- theta[2]
c <- theta[3]
tau.sq <- theta[4]
rnorm(length(x), fn(a,b,c,x), sqrt(tau.sq))
}
post.curves <- t(apply(samples.1[burn.in:nrow(samples.1),], 1, fn.pred, x))
post.curves.quants <- summary(mcmc(post.curves))$quantiles
plot(x, y, pch=19, xlab="x", ylab="f(x)")
lines(x, post.curves.quants[,1], lty="dashed", col="blue")
lines(x, post.curves.quants[,3])
lines(x, post.curves.quants[,5], lty="dashed", col="blue")
## End(Not run)[＃（不执行）]
转载请注明:出自 生物统计家园网(http://www.biostatistic.net)。
注：
注1：为了方便大家学习，本文档为生物统计家园网机器人LoveR翻译而成，仅供个人R语言学习参考使用，生物统计家园保留版权。
注2：由于是机器人自动翻译，难免有不准确之处，使用时仔细对照中、英文内容进行反复理解，可以帮助R语言的学习。
注3：如遇到不准确之处，请在本贴的后面进行回帖，我们会逐渐进行修订。
|
发表于 2012-9-30 14:22:18
