1 #include "bayesm.h"
2
3 //[[Rcpp::export]]
breg(vec const & y,mat const & X,vec const & betabar,mat const & A)4 vec breg(vec const& y, mat const& X, vec const& betabar, mat const& A) {
5
6 // Keunwoo Kim 06/20/2014
7
8 // Purpose: draw from posterior for linear regression, sigmasq=1.0
9
10 // Output: draw from posterior
11
12 // Model: y = Xbeta + e e ~ N(0,I)
13
14 // Prior: beta ~ N(betabar,A^-1)
15
16 int k = betabar.size();
17 mat RA = chol(A);
18 mat W = join_cols(X, RA); //same as rbind(X,RA)
19 vec z = join_cols(y, RA*betabar);
20 mat IR = solve(trimatu(chol(trans(W)*W)), eye(k,k)); //trimatu interprets the matrix as upper triangular and makes solve more efficient
21
22 return ((IR*trans(IR))*(trans(W)*z) + IR*vec(rnorm(k)));
23 }
24