请问使用metrop要注意什么 谢谢

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

请问除了function的形式必须是 log unnormalized probability density 以及初值代入后大于零 还有什么要注意的吗？
log unnormalized probability density 是不是就是把原有的function外面加个log()就好了？谢谢~~~

我的始终报错如下：
错误于system.time(out <- .Call("metrop", func1, initial, nbatch, blen,  :
  logh: func returned NA or NaN

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏0 回帖

关键词：Metro Etro Probability bability function function

metrop(obj, initial, nbatch, blen = 1,.....)
  obj : an R function that evaluates the log unnormalized probability density
        of the desired equilibrium distribution of the Markov chain.
  return value 可以是 -Inf
  但不可以是Inf, NA or NaN

所以底下范例是OK
  h <- function(x) if (all(x >= 0) && sum(x) <= 1) return(1) else return(-Inf)
  out <- metrop(h, rep(0, 5), 1000)

若改成底下就出错了
  h <- function(x) if (all(x >= 0) && sum(x) <= 1) return(1) else return(Inf)
  out <- metrop(h, rep(0, 5), 1000)

epoh 发表于 2012-4-20 21:39
metrop(obj, initial, nbatch, blen = 1,.....)
  obj : an R function that evaluates the log unnormali ...

谢谢~~请问所以是不是我还需要再加个判断语句判断h的可能结果?当我不用metrop 把向量i赋值代入方程h是有一个值的~
try<-c(a,b,c)
h<-function(try){
atoxity<-0.05
btoxity<-0.1
n<-1
x<-0
pie<-1-((1-atoxity^try[1])^(-try[3])+(1-btoxity^try[2])^(-try[3])-1)^(-1/try[3])
logl<-log((pie^x)*(1-pie)^(n-x))
logprior<-log(dgamma(try[1],0.5,0.5)*dgamma(try[2],0.5,0.5)*dgamma(try[3],0.1,0.1))
return (logl+logprior)
}
library(mcmc)
set.seed(528)
i<-c(3.334313e-01,5.218204e-02,2.324598e-14)
out <- metrop(h,i,100,blen = 1, nspac = 1, scale = 0.5)

ccjifei 发表于 2012-4-21 00:09
谢谢~~请问所以是不是我还需要再加个判断语句判断h的可能结果?当我不用metrop 把向量i赋值代入方程h是有一 ...

你可以参考MCMC Package Example
demo.pdf来写log unnormalized posterior (log likelihood plus log prior) density function
#####
library(mcmc)
data(logit)
out <- glm(y ~ x1 + x2 + x3 + x4, data = logit, family = binomial(), x = TRUE)
summary(out)
#log unnormalized posterior (log likelihood plus log prior) density
lupost <- function(beta, x, y) {
eta <- as.numeric(x %*% beta)
logp <- ifelse(eta < 0, eta - log1p(exp(eta)), -log1p(exp(-eta)))
logq <- ifelse(eta < 0, -log1p(exp(eta)), -eta - log1p(exp(-eta)))
logl <- sum(logp[y == 1]) + sum(logq[y == 0])
return(logl - sum(beta^2)/8)
}

set.seed(42)
x <- out$x
y <- out$y
beta.init <- as.numeric(coefficients(out))
metropout <- metrop(lupost, beta.init, 1000, x = x, y = y)
names(metropout)

谢谢啊 实在不好意思 我看了demo 想请教一下 这两句话是什么意思呀 log unnormalized posterior density function 是不是就是函数外加上Log然后写成log likelihood plus log prior的样子呢 真心麻烦你了~~
logp <- ifelse(eta < 0, eta - log1p(exp(eta)), -log1p(exp(-eta)))
logq <- ifelse(eta < 0, -log1p(exp(eta)), -eta - log1p(exp(-eta)))

ccjifei 发表于 2012-4-22 00:38
谢谢啊 实在不好意思 我看了demo 想请教一下 这两句话是什么意思呀 log unnormalized posterior density fu ...

#------------------
# data
#------------------

y <- c(-0.91, 1.27, 0.85, 0.89, -0.22, 1.07, 2.25, -1.20, -0.95, -0.82)
n <- length(y)

#-----------------------------------------------------------------
# model states that
#
# [Y | mu, psi] i.i.d~ Normal(mu, 1 / psi),   i = 1, ..., n
#
# Parameter vector theta = (mu, psi), where mu is the expected value
# and psi is the precision of the normal population.
#
# use the improper prior with density proportional to 1/psi.
#--------------------------------------------------------------------
# Define functions for calculating  log-likelihood, log-prior, and
# log(unnormalized posterior)
#--------------------------------------------------------------------

# Log-likelihood
loglik <- function(theta) {
  sum(dnorm(log = TRUE, y, mean = theta[1], sd = 1 / sqrt(theta[2])))
}

# log-prior for the improper prior
logprior <- function(theta) (-log(theta[2]))

# Log unnormalized posterior
logpost <- function(theta) (loglik(theta) + logprior(theta))


##############
ifelse用法请看帮助
ifelse(test, yes, no)
     test an object which can be coerced to logical mode.
     yes return values for true elements of test.
     no return values for false elements of test.

logp <- ifelse(eta < 0, eta - log1p(exp(eta)), -log1p(exp(-eta)))

当 eta < 0(yes),logp = eta - log1p(exp(eta))
否则      (no) ,logp =  -log1p(exp(-eta))