机器算法验证 - 在 R 中实现 Metropolis Hastings 算法 - 吾爱随笔录

在 R 中实现 Metropolis Hastings 算法

机器算法验证自习贝叶斯算法蒙特卡洛大都会黑斯廷斯

2022-04-09 12:03:03

考虑一个具有均值和方差的单变量正态模型。使用 Beta(2,2) 先验（不知何故，我们知道 µ 在 0 和 1 之间），对先验（回想一下，如果随机变量是那么是。先验假设和是独立的。实施 Metropolis-Hastings 算法来评估和的后验分布。请记住，您必须共同接受或拒绝和 $µ$ $τ$ $µ$ $log-normal(1,10)$ $τ$ $X$ $log-normal(m,v)$ $log X$ $N(m,v))$ $µ$ $τ$ $µ$ $τ$ $µ$ $τ$ . $µ$ 大于 0.5的后验概率。

以下是数据：

2.3656491 2.4952035 1.0837817 0.7586751 0.8780483 1.2765341
1.4598699 0.1801679 -1.0093589 1.4870201 -0.1193149 0.2578262

尝试：（这里输入我的数据）

mu <- rbeta(1,2,2)
tau <- rlnorm(1,1,10)
x <- 
  c(2.3656491, 2.4952035, 1.0837817, 0.7586751, 0.8780483, 1.2765341,
  1.4598699, 0.1801679, -1.0093589, 1.4870201, -0.1193149, 0.2578262)

这是网上找到的MH算法。我不确定如何将上述数据应用于它。

# starting point
a0 <- -5 
b0 <- -10 

# length of chain
nit <- 1000
a <- rep(0,nit)
b <- rep(0,nit)

# initialize
a[1] <- a0
b[1] <- b0

# tuning parameter
s0 <- 2.0 # maximum step size in random walk proposal
# function. try different s0, e.g., 0.1, 1.0, 2.0

# start chain
counter = 0 # monitor number of acceptances.
for( i in 1:nit) {
  s <- s0*runif(1)
  theta<-2*3.1415926*runif(1)
  anew <- a[i] + s*cos(theta) # random walk
  bnew <- b[i] + 2*s*sin(theta) # random walk 
  r <- pdf(anew,bnew)/pdf(a[i],b[i])# acceptance ratio
  test <- runif(1)
  if(test < r ) # accept proposed moved.
  {     
     a[i+1] <- anew
     b[i+1] <- bnew
    counter = counter + 1;
  }
  else # reject proposed move, stay put.
  {
    a[i+1] <- a[i]
    b[i+1] <- b[i]
  }
}
# acceptance rate =
counter /nit

1个回答

仅仅从没有解释的在线代码开始学习一个主题确实是一个非常糟糕的主意。最好阅读一本书（比如我们的蒙特卡洛方法与 R 简介！）或介绍性论文并编写自己的代码。

进行随机游走，在您的情况下可能是以说明与随机游走提议不兼容和这个随机游走建议是从以下行推导出来的 $(a,b)$

(logit (μ), \log (τ)) = (\log {\frac{μ}{1 - μ}}, \log (τ))

$(\text{logit}(\mu),\log(\tau))=\left(\log\left\{\dfrac{\mu}{1-\mu}\right\},\log(\tau)\right)$

μ

$\mu$

τ

$\tau$

  theta<-2*3.1415926*runif(1)
  anew <- a[i] + s*cos(theta) # random walk
  bnew <- b[i] + 2*s*sin(theta) # random walk

这对应于模拟二元法线（直接调用也可以！）并将第一个组件缩放和第二个组件。（没有理由造成这种差异，每个组件都需要基于接受率的自己的规模。）这条线 $\text{N}(0,1)$ rnorm(2)s2s

s <- s0*runif(1)

在每一步改变比例实际上是不正确的，因为它使通常的 Metropolis-Hasting 比率无效：在这个比率中

r <- pdf(anew,bnew)/pdf(a[i],b[i])# acceptance ratio

pdf表示目标分布的密度，在您的情况下，可以编码为

pdf <-function(a,b){

  dbeta(1/(1+exp(-a)),2,2)*dnorm(b,1,10)*prod(
  dnorm(x,mean=1/(1+exp(-a)),sd=exp(b)))*
  exp(a)/(1+exp(a))^2}

它将先验乘以雅可比行列式的可能性。这听起来可能很复杂，但解释是

随机游走对参数化进行操作，Metropolis-Hastings 比率是该参数化目标的比率； $(\text{logit}(\mu),\log(\tau))$
先验是为参数化或的，并且必须转换为相同参数化的先验上使用 Normal 模型，上使用 Beta(2,2) 先验的 logit 变换，因此逆 logit 变换的雅可比： $(\mu,\tau)$ $(\mu,\log\,\tau)$ $(\text{logit}(\mu),\log(\tau))$ $\log(\tau))$ $\mu$ $\frac{\exp (a)}{(1 + \exp (a))^{2}}$ $\dfrac{\exp(a)}{(1+\exp(a))^2}$

一旦你定义了这个函数，你的程序应该可以运行，至少从理论上来说是这样。如果接受率太高或太低，则需要校准秤s。

其它你可能感兴趣的问题

上一篇Hyndman 对适当时间序列交叉验证的解释与 Leave-One-Out 有何不同？下一篇LIBSVM“警告：使用 -h 0 可能更快”