代码之家 › 专栏 › 技术社区 › EDennnis

使用Apply族的迭代算法:Metropolis-Hasting-Alg

lapply for-loop algorithm r

EDennnis · 技术社区 · 7 年前

#For Loop Application of Metropolis-Hastings
phi <- matrix(,m+1,1) # (m+1) x 1 vector to save samples of x in
phi[1,] <- 0.5 # Starting value

set.seed(1603)
for(k in 1:m){
 phi.star <- runif(1) # Uniform proposal distribution
 phi.k <- phi[k,1]
 R <- min(1,(dbinom(1,1,phi.star)*dbeta(phi.star,1,1))/(dbinom(1,1,phi.k)*dbeta(phi.k,1,1)))
 phi[k+1,] <- dplyr::case_when(R > runif(1) ~ phi.star,
                            TRUE ~ phi.k) #Retain phi.star with probability R
}

#Lapply Function for Metropolis-Hastings
phi <- list()
set.seed(1603)
phi <- lapply(X = 1:m, function(k){
  phi.star <- runif(1) # Uniform proposal distribution
  phi.k <- phi[[k]]
  R <- min(1,(dbinom(1,1,phi.star)*dbeta(phi.star,1,1))/(dbinom(1,1,phi.k)*dbeta(phi.k,1,1)))
  phi.keep <- dplyr::case_when(R > runif(1) ~ phi.star,
                           TRUE ~ phi.k) #Retain phi.star with probability R
})

phi <- phi %>% do.call("rbind",.)

出现这个问题是因为我试图访问上一个列表值。任何帮助都将不胜感激。

1 回复 | 直到 7 年前

duckmayr 7 年前

下面是我的评论的一个答案版本,以及如何实现它的一个示例,并对性能提高进行了比较。

apply() 不是路要走; you really need to loop . 如果您想加速循环,一个好的方法是转到编译代码,这可以通过 Rcpp . 从哈德利那里 Advanced R “C++所能解决的典型瓶颈包括:循环不易被量化,因为后续迭代依赖于以前的迭代……”

在演示如何在C++中实现 Rcpp公司 phi.star/phi.k 自从 dbinom(1, 1, x) = x 十 dbeta(x, 1, 1) =1表示任何 十 .

以下是在Rcpp中实现M-H采样器的一种方法:

#include <Rcpp.h>

using namespace Rcpp;

// [[Rcpp::export]]
NumericVector mh_cpp(double starting_value, int n) {
    NumericVector phi(n+1);
    phi[0] = starting_value;
    for ( int i = 0; i < n; ++i ) {
        double phi_star = R::runif(0.0, 1.0);
        double phi_k = phi[i];
        double r = phi_star / phi_k;
        if ( r >= 1 || R::runif(0.0, 1.0) < r ) {
            phi[i+1] = phi_star;
        }
        else {
            phi[i+1] = phi_k;
        }
    }
    return phi;
}

然后我为你的采样器创建了一个R函数(注意代码有多相似!),并比较了性能:

mh_r <- function(starting_value, n) {
    phi <- numeric(n+1)
    phi[1] <- starting_value
    for ( k in 1:n ) {
        phi_star <- runif(1) # Uniform proposal distribution
        phi_k <- phi[k]
        r <- phi_star / phi_k
        phi[k+1] <- "if"(r >= 1 | runif(1) < r, phi_star, phi_k)
    }
    return(phi)
}

m <- 100000
library(microbenchmark)
microbenchmark(mh_r(0.5, m), mh_cpp(0.5, m))

Unit: milliseconds
           expr        min         lq      mean     median         uq        max
   mh_r(0.5, m) 2355.68150 2376.59047 2421.6640 2383.96823 2408.27571 3816.37139
 mh_cpp(0.5, m)   10.54044   10.59464   10.8235   10.61732   10.65326   25.43983

m <- 1000
microbenchmark(mh_r(0.5, m), mh_cpp(0.5, m))
Unit: microseconds
           expr       min         lq       mean     median        uq       max
   mh_r(0.5, m) 22199.272 22395.6200 24223.8546 22511.8160 22705.792 39525.690
 mh_cpp(0.5, m)   115.186   118.4795   130.4884   131.0385   135.016   403.093

Rcpp公司 ,在1000次迭代和100000次迭代中,我们得到了大约200倍的速度。

Rcpp公司 ,我建议 Rcpp Introduction Vignette 哈德利的那一章高级R