Description Usage Arguments Details Value Examples
Demonstrates the use of RcppDist in C++ with Bayesian linear regression
1 | bayeslm(y, x, iters = 1000L)
|
y |
A numeric vector – the response |
x |
A numeric matrix – the explanatory variables; note this assumes you have included a column of ones if you intend there to be an intercept. |
iters |
An integer vector of length one, the number of posterior draws desired; the default is 1000. |
To see an example of using RcppDist C++ functions in C++ code, we can code up a Bayesian linear regression with completely uninformative priors (such that estimates should be equivalent to classical estimates). The code to do so is as follows:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 | #include <RcppDist.h>
// or, alternatively,
// #include <RcppArmadillo.h>
// #include <mvnorm.h>
// [[Rcpp::depends(RcppArmadillo, RcppDist)]]
// [[Rcpp::export]]
Rcpp::List bayeslm(const arma::vec& y, const arma::mat x,
const int iters = 1000) {
int n = x.n_rows;
int p = x.n_cols;
double a = (n - p) / 2.0;
arma::mat xtx = x.t() * x;
arma::mat xtxinv = xtx.i();
arma::vec mu = xtxinv * x.t() * y;
arma::mat px = x * xtxinv * x.t();
double ssq = arma::as_scalar(y.t() * (arma::eye(n, n) - px) * y);
ssq *= (1.0 / (n - p));
double b = 1.0 / (a * ssq);
arma::mat beta_draws(iters, p);
Rcpp::NumericVector sigma_draws(iters);
for ( int iter = 0; iter < iters; ++iter ) {
double sigmasq = 1.0 / R::rgamma(a, b);
sigma_draws[iter] = sigmasq;
// Here we can use our multivariate normal generator
beta_draws.row(iter) = rmvnorm(1, mu, xtxinv * sigmasq);
}
return Rcpp::List::create(Rcpp::_["beta_draws"] = beta_draws,
Rcpp::_["sigma_draws"] = sigma_draws);
}
|
A list of length two; the first element is a numeric matrix of the beta draws and the second element is a numeric vector of the sigma draws
1 2 3 4 5 6 7 8 9 10 |
[1] 9.89 2.05 -0.91 3.01
[1] 9.89 2.04 -0.91 3.02
[1] 10 2 -1 3
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.