import ramcmc

DESCRIPTION

@@ -16,7 +16,8 @@ LinkingTo:
     bsvars, 
     Rcpp, 
     RcppArmadillo, 
-    RcppProgress
+    RcppProgress,
+    ramcmc
 Depends: 
     R (>= 2.10),
     bsvars

src/mcmc.cpp

@@ -0,0 +1,70 @@
+
+#include <functional>
+#include <iostream>
+#include <RcppArmadillo.h>
+#include "Rcpp/Rmath.h"
+#include "progress.hpp"
+#include "ramcmc.h"
+
+using namespace arma;
+
+
+// adaptive Metropolis algorithm for strictly positive parameters
+// [[Rcpp:interfaces(cpp)]]
+arma::mat metropolis(
+    const int& T,
+    const int& t0,
+    arma::vec  x,
+    arma::mat  Sigma,
+    const std::function<double(const arma::vec&)>& log_target
+) {
+
+  int    n = x.n_elem;
+
+  x        = log(x);
+  vec xbar = x;
+  double s = 2.38 / sqrt(n);
+  double d = log_target(x) + sum(x);
+
+  double new_d, a;
+  vec    new_x, diff;
+
+  mat X(n, T);
+  X.col(0) = x;
+
+  vec      progress = round(linspace(0, T, 50));
+  Progress p(50, true);
+
+  for (int t = 1; t < T; t++) {
+    new_x = mvnrnd(x, s*s * Sigma);
+    new_d = log_target(exp(new_x)) + sum(new_x);
+    a     = std::min(1.0, exp(new_d - d));
+
+    if (randu() < a) {
+      x = new_x;
+      d = new_d;
+    }
+    X.col(t) = x;
+
+    if (t <= t0) {
+      diff   = x - xbar;
+      s     += pow(t, -0.6) * (a - 0.234);
+      xbar  += diff / (t + 1);
+      Sigma  = Sigma * t / (t + 1) + diff * diff.t() * t / (t + 1) / (t + 1);
+    }
+
+    if (any(progress == t)) {
+      p.increment();
+    }
+  }
+
+  X = exp(X);
+  return X;
+}
+
+
+
+
+
+
+

src/mcmc.h

@@ -0,0 +1,23 @@
+#ifndef _MCMC_H_
+#define _MCMC_H_
+
+
+#include <functional>
+#include <iostream>
+#include <RcppArmadillo.h>
+#include "Rcpp/Rmath.h"
+#include "progress.hpp"
+
+using namespace arma;
+
+
+arma::mat metropolis(
+    const int& T,
+    const int& t0,
+    arma::vec  x,
+    arma::mat  Sigma,
+    const std::function<double(const arma::vec&)>& log_target
+);
+
+
+#endif  // _MCMC_H_

src/sample_hyper.cpp

@@ -2,6 +2,7 @@
 #include <RcppArmadillo.h>
 
 #include "utils.h"
+#include "mcmc.h"
 
 using namespace Rcpp;
 using namespace arma;

src/utils.cpp

@@ -2,8 +2,6 @@
 #include <functional>
 #include <iostream>
 #include <RcppArmadillo.h>
-#include "Rcpp/Rmath.h"
-#include "progress.hpp"
 
 using namespace arma;
 
@@ -79,66 +77,3 @@ arma::mat Df(
   return result;
 }
 
-
-// adaptive Metropolis algorithm for strictly positive parameters
-// [[Rcpp:interfaces(cpp)]]
-arma::mat metropolis(
-    const int& T,
-    const int& t0,
-    arma::vec  x,
-    arma::mat  Sigma,
-    const std::function<double(const arma::vec&)>& log_target
-) {
-
-  int    n = x.n_elem;
-
-         x = log(x);
-  double s = 2.38 / sqrt(n);
-  double d = log_target(x) + sum(x);
-
-  double new_d, a;
-  vec    new_x, xbar, diff;
-
-  mat X(n, T);
-  X.col(0) = x;
-
-  vec      progress = round(linspace(0, T, 50));
-  Progress p(50, true);
-
-  for (int t = 1; t < T; t++) {
-    new_x = mvnrnd(x, s*s * Sigma);
-    new_d = log_target(exp(new_x)) + sum(new_x);
-    a     = std::min(1.0, exp(new_d - d));
-
-    if (randu() < a) {
-      x = new_x;
-      d = new_d;
-    }
-    X.col(t) = x;
-
-
-    if (t == 1) {
-      xbar   = mean(X.cols(0, t), 1);
-      Sigma  = cov(X.cols(0, t).t());
-    } else if (t <= t0) {
-      diff   = x - xbar;
-      s     += pow(t, -0.6) * (a - 0.234);
-      xbar  += diff / (t + 1);
-      Sigma  = Sigma * t / (t + 1) + diff * diff.t() * t / (t + 1) / (t + 1);
-    }
-
-    if (any(progress == t)) {
-      p.increment();
-    }
-  }
-
-  X = exp(X);
-  return X;
-}
-
-
-
-
-
-
-