Skip to content

Baby steps, continued

March 19, 2013

I’ve found some more bugs in my code without even looking for them:

I added an entry point for simulating from the Potts model for a fixed value of β, without any external field. This function can be used to estimate the critical value of β by simulation (since β doesn’t depend directly on the data). In doing this, I discovered that my neighbourhood matrix wasn’t a matrix … it was a vector! It appears that the following code silently converts from Rcpp::IntegerMatrix to Rcpp::IntegerVector and thence into an arma::umat with only a single column:

arma::umat neigh = unsign(nR - 1);

This is because Rcpp::Matrix inherits from Rcpp::Vector (shouldn’t it be the other way around, generalization not specialization?) and the binary arithmetic operators are defined to return Rcpp::Vector. Exactly the kind of shenanigans that Felix von Leitner’s polemic warns about…

So, once I fixed bug number 1, this uncovered other problems deeper inside my code:

double total_prob = 0.0;
const arma::uvec neighbors = neigh.row(i);
const arma::umat labels = z.rows(neighbors);
for (unsigned j=0; j < labels.n_cols; j++)
   unsigned sum_neigh = arma::sum(labels.col(j));
   total_prob += exp(beta*sum_neigh);

The  above code doesn’t even come close to performing the vectorized operations that I was expecting. In the end, I gave up and just wrote it as a nested loop:

double total_prob = 0.0;
for (unsigned j=0; j < z.n_cols; j++)
   unsigned sum_neigh = 0;
   for (unsigned k=0; k < neigh.n_cols; k++)
      sum_neigh += z(neigh(i,k),j);
   total_prob += exp(beta*sum_neigh);

Trying to be too clever just hasn’t paid off with RcppArmadillo.

From → C++, MCMC

One Comment

Trackbacks & Pingbacks

  1. Critical Point | Matt Moores

Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Ella Kaye on Ella Kaye

Computational Bayesian statistics

Bayes' Food Cake

A bit of statistics, a bit of cakes. - Blogs to Learn R from the Community

Computational Bayesian statistics

Richard Everitt's blog

Computational Bayesian statistics

Let's Look at the Figures

David Firth's blog

Nicholas Tierney

Computational Bayesian statistics

Mad (Data) Scientist

Musings, useful code etc. on R and data science

Another Astrostatistics Blog

The random musings of a reformed astronomer ...

Darren Wilkinson's research blog

Statistics, computing, data science, Bayes, stochastic modelling, systems biology and bioinformatics

(badness 10000)

Computational Bayesian statistics

Igor Kromin

Computational Bayesian statistics


I can't get no

Xi'an's Og

an attempt at bloggin, nothing more...

Sam Clifford

Postdoctoral Fellow, Bayesian Statistics, Aerosol Science

%d bloggers like this: