[Rcpp-devel] CppBugs vs WinBugs

Watson, Samuel S.I.Watson at warwick.ac.uk
Mon Dec 19 14:10:05 CET 2011

```The radon model is as follows (I haven’t changed the linear model posted on Github), I slightly altered the radon model from Github which is marked inline below:

BUGS code (to be saved in WinBugs working directory as "radon.bug"):

model {
for (i in 1:919){
y.hat[i] <- a[county2[i]] + b*floor[i]
}
b ~ dnorm (0, .0001)
tau.y <- pow(sigma.y, -2)
sigma.y ~ dunif (0, 100)

for (j in 1:85){
a[j] ~ dnorm (mu.a, tau.a)
}
mu.a ~ dnorm (0, .0001)
tau.a <- pow(sigma.a, -2)
sigma.a ~ dunif (0, 100)
}
-------------------------------------
R code:

require(inline)
require(Rcpp)
require(R2WinBUGS)

model<-'
#include <iostream>
#include <fstream>
#include <vector>
#include <string>
#include <algorithm>
#include <cmath>
#include <boost/random.hpp>
#include <boost/algorithm/string.hpp>
#include <cppbugs/cppbugs.hpp>

using namespace arma;
using namespace cppbugs;
using std::vector;
using std::string;
using std::cout;
using std::endl;
using std::ifstream;

public:
const vec& level;
const vec& basement;
const mat& group;
int N, N_counties;
mat indicator_matrix; //indicator matrix used for county level random effects

Normal<vec> a;
Normal<double> b;
Deterministic<double> tau_y;
Uniform<double> sigma_y;
Normal<double> mu_a;
Deterministic<double> tau_a;
Uniform<double> sigma_a;
Deterministic<mat> y_hat;
Normal<mat> likelihood;

RadonVaryingInterceptModel(const vec& level_, const vec& basement_, const mat& group_,int N_, int N_counties_):
level(level_),basement(basement_),group(group_),
a(randn<vec>(N_counties_)), b(0),N(N_),N_counties(N_counties_),
indicator_matrix(N,N_counties),
tau_y(1),sigma_y(1),mu_a(0),tau_a(1),sigma_a(1),
y_hat(randn<mat>(level_.n_rows,1)),likelihood(level_,true)

{
indicator_matrix.fill(0.0);
for(int i=0;i<group.n_elem;i++){
indicator_matrix(i,group[i])=1.0;
}
}

void update() {
y_hat.value = indicator_matrix*a.value + b.value*basement;
tau_y.value = pow(sigma_y.value, -2.0);
tau_a.value = pow(sigma_a.value, -2.0);
a.dnorm(mu_a.value, tau_a.value);
b.dnorm(0, 0.0001);
sigma_y.dunif(0, 100);
mu_a.dnorm(0, 0.0001);
sigma_a.dunif(0, 100);
likelihood.dnorm(y_hat.value,tau_y.value);
}
};
'
Src<-'

mat group = Rcpp::as<arma::mat>(GR);
vec level = Rcpp::as<arma::vec>(lev);
vec basement = Rcpp::as<arma::vec>(basm);
int N = 919;
int N_counties = 85;
int iterations_ = as<int>(iterations);
int burn_ = as<int>(burn);
int thin_ = as<int>(thin);

return Rcpp::List::create(Rcpp::Named("b", m.b.mean()), Rcpp::Named("ar", m.acceptance_ratio()),
Rcpp::Named("a", m.a.mean()));
'

df.mn <- df[df[,2]=="MN",]
floor <- df.mn\$floor

# get county index variable
county.name <- as.vector(df.mn\$county)
uniq <- unique(county.name)
J <- length(uniq)
county <- rep (NA, J)
for (i in 1:J){
county[county.name==uniq[i]] <- i
}

#its necessary to change the structure of the data files or else winbugs has an error
floor<-as.numeric(floor)
county<-as.numeric(county)

#n.burnin = n.iter/2 for bugs

-----------------------------------------------------------------------
The difference I get is 16.53s vs 13.85s

Kind regards,
Sam Watson

-----Original Message-----
From: Whit Armstrong [mailto:armstrong.whit at gmail.com]
Sent: 19 December 2011 12:39
To: Watson, Samuel
Cc: rcpp-devel at r-forge.wu-wien.ac.at
Subject: Re: [Rcpp-devel] CppBugs vs WinBugs

The speedup depends on  your model.  My comparisons were based on R/JAGS and pymc on linux.

-Whit

On Mon, Dec 19, 2011 at 7:30 AM, Watson, Samuel <S.I.Watson at warwick.ac.uk> wrote:
> I am currently using CppBugs with Rcpp through R. I am very interested
> to use CppBugs as I am finding WinBugs to be prohibitively slow, I use
> large amounts of data in large multilevel models, so when I found
> cppbugs I was excited. It says on the Github page for cppbugs that I
> can achieve speeds of 20-100x faster than winbugs.
>
> I have been testing some examples: for the linear model example that
> is provided at Github, the time difference is about 2x
>
> I also tested the Radon model with Rcpp and Andrew Gelmans data and
> for
> 10000 iterations and 10000 burnin with a thinning parameter of 5 the
> difference is only 16 seconds vs 13 seconds.
>
>
>
> I can run multiple chains in parallel through R with the ‘snow’
> package,
> cxxfunction() and package.skeleton() and this does seem to provide a
> little bit of a boost (parallel winbugs vs parallel cppbugs). But
> nothing greater than 4x.
>
>
>
> Is it possible to achieve the kind of difference mentioned on Whit
> Armstrong’s page for cppbugs with R?
>
>
>
> I am happy to post all the code if required.
>
>
>
> Many thanks
>
> Sam Watson
>
>
> _______________________________________________
> Rcpp-devel mailing list
> Rcpp-devel at lists.r-forge.r-project.org
> https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-deve
> l
```