r - 从概率网格中采样（贝叶斯后验近似）答案

【问题标题】：r - Sampling from a grid of probabilities (Bayesian posterior approximation)r - 从概率网格中采样（贝叶斯后验近似）
【发布时间】：2015-05-03 02:56:51
【问题描述】：

我正在做贝叶斯分析，并试图估计两个参数。为了近似后验分布，我构建了一个精细网格并计算了网格中每个元素的后验概率。我对其进行了标准化，使网格总和为 1。

现在我对从分布中抽样感兴趣。这是我目前所拥有的：

sampleGrid <- function(post.grid, mu.grid, sig2.grid) {
  value <- sample(post.grid, 1, prob=post.grid)
  index <- which(post.grid == value) 
  col <- as.integer(index/nrow(post.grid))+1
  row <- index-(col-1)*nrow(post.grid)
  return(c(mu.grid[row], sig2.grid[col]))
}

但是，当我想进行大量采样时，我遇到了运行时问题，因为我使用了 for 循环：

for(i in 1:nrow(sample.grid)) {
  sample.grid[i, ] <- sampleFromGrid(post.grid, mu.grid, sig2.grid)
}

我想知道是否有办法对此进行矢量化。我的尝试是：

vectorizedSampleFromGrid <- function(post.grid, mu.grid, sig2.grid, n){
    values <- sample(post.grid, n, replace=T, prob=post.grid)
    index <- which(post.grid %in% values)
    if(length(values)!=length(index)) {
        temp.df <- count(values)
        index <- which(post.grid %in% temp.df[,1])
        temp.df <- cbind(temp.df, index)
        temp.df <- temp.df[temp.df[, 2] > 1, ]
        for(i in 1:nrow(temp.df)) {
            index <- c(index, rep(temp.df[i, 3], temp.df[i,2]-1))
        }
    }
    col <- as.integer(index/nrow(post.grid))+1
    row <- index-(col-1)*nrow(post.grid)
    return(cbind(mu.grid[row], sig2.grid[col]))
}

我知道有些元素会被多次采样。我想要做的是根据它们被采样的次数将这些索引多次附加到原始索引列表中。但是，当我这样做时，结果不正确。

如果有人能提供任何建议，我将不胜感激。

【问题讨论】：

标签： r bayesian random-sample approximation

【解决方案1】：

这就是我要做的。创建一个矢量化函数来评估后验（或至少与其成比例的东西）：

f = function(mu, sigma, log=TRUE) {
  logf = dnorm(mu, 0, sigma, log=TRUE) + dgamma(sigma, 1, 1, log=TRUE)
  if (log) return(logf)
  return(exp(f))
}

现在在网格上评估这个函数。

library(dplyr)
grid = mutate(expand.grid(mu=seq(-3,3,1), sigma=seq(1,7,1)),
              logp = f(mu,sigma),
              logp = logp-max(logp), # for numerical stability
              p    = exp(logp),
              p    = p/sum(p))       # Normalize

现在从这个网格中获取样本：

samples = sample_n(grid, size=100, replace=TRUE, weight=grid$p)

【讨论】：