r - 從概率網格抽樣（貝葉斯後驗近似）

我正在做貝葉斯分析，我試圖估計兩個參數。爲了近似後驗分佈，我構建了一個精細網格並計算網格中每個元素的後驗概率。我規範化它，使網格總和爲1.r - 從概率網格抽樣（貝葉斯後驗近似）

現在我對分佈採樣感興趣。這是我到目前爲止有：

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])) 
}

我知道一些元素將被採樣一次以上。我試圖做的是將這些索引多次添加到原始索引列表中，取決於它們被採樣了多少次。但是，當我這樣做時，結果是不正確的。

如果有人能提供任何建議，我將不勝感激。

來源

2015-03-02 user3821273

這是我會做的。創建一個矢量化函數來評估後驗（或至少與其成正比的東西）：

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)

來源

2015-03-04 16:03:39 jaradniemi

r - 從概率網格抽樣（貝葉斯後驗近似）

回答

相關問題