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.

Advertisements

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:

WordPress.com Logo

You are commenting using your WordPress.com 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

Let's Look at the Figures

David Firth's blog

Nicholas Tierney

Computational Bayesian statistics

One weiRd tip

Computational Bayesian statistics

Series B'log

discussion blog for JRSS Series B papers

Mad (Data) Scientist

Musings, useful code etc. on R and data science

R-bloggers

R news and tutorials contributed by (750) R bloggers

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

CHANCE

Computational Bayesian statistics

StatsLife - Significance magazine

Computational Bayesian statistics

(badness 10000)

Computational Bayesian statistics

Igor's Blog

Computational Bayesian statistics

Statisfaction

I can't get no

Xi'an's Og

an attempt at bloggin, nothing more...

Sam Clifford

Postdoctoral Fellow, Bayesian Statistics, Aerosol Science

Bayesian Research & Applications Group

Frontier Research in Bayesian Methodology & Computation

%d bloggers like this: