Итерационные алгоритмы с использованием семейства Apply: Метрополис-Хастинг Альг

Question

При построении алгоритма нам часто нужны как предыдущие, так и текущие значения итерации.Я пытаюсь использовать семейство применений, чтобы улучшить скорость алгоритма.Я построил алгоритм Метрополиса-Хастинга с использованием цикла for, но мне нужна помощь с использованием семейства apply (я предпочитаю использовать lapply).Я приложил код, который я построил до сих пор.

#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",.)

Проблема возникает из-за того, что я пытаюсь получить доступ к предыдущему значению списка.Любая помощь будет оценена.

Answer 1

Вот вариант ответа на мой комментарий с примером того, как его реализовать, и сравнение увеличения производительности.

Поскольку вы полагаетесь на предыдущее значение, apply() непуть; вам действительно нужно цикл .Если вы хотите ускорить циклы, хороший способ - перейти к скомпилированному коду, что довольно легко сделать с помощью Rcpp.Из Advanced R от Хэдли: «Типичные узкие места, которые может устранить C ++, включают в себя: Циклы, которые нельзя легко векторизовать, поскольку последующие итерации зависят от предыдущих…» *

Прежде чем показать, какреализуя в C ++ через Rcpp, я сначала отмечу, что ваше отношение MH здесь упрощается до phi.star/phi.k, поскольку dbinom(1, 1, x) = x для любых x и dbeta(x, 1, 1) = 1 для любых x.

Вот один из способов реализации вашего сэмплера MH в Rcpp:

#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

Итак, перейдя на C ++ с помощью Rcpp, мы получили, что это будет примерно в 200 раз быстрее как при 1000 итерациях, так и при 100000 итерациях.

Для некоторой вводной информации об использовании Rcpp, Я бы предложил Rcpp Введение Виньетка и главу Хэдли об этом в Advanced R .