C++ libraries for statistical computing

Question

I've got a particular MCMC algorithm which I would like to port to C/C++. Much of the expensive computation is in C already via Cython, but I want to have the whole sampler written in a compiled language so that I can just write wrappers for Python/R/Matlab/whatever.

After poking around I'm leaning towards C++. A couple of relevant libraries I know of are Armadillo (http://arma.sourceforge.net/) and Scythe (http://scythe.wustl.edu/). Both try to emulate some aspects of R/Matlab to ease the learning curve, which I like a lot. Scythe squares a little better with what I want to do I think. In particular, its RNG includes a lot of distributions where Armadillo only has uniform/normal, which is inconvenient. Armadillo seems to be under pretty active development while Scythe saw its last release in 2007.

So what I'm wondering is if anyone has experience with these libraries -- or others I have almost surely missed -- and if so, whether there is anything to recommend one over the others for a statistician very familiar with Python/R/Matlab but less so with compiled languages (not completely ignorant, but not exactly proficient...).

Answer 1

We have spent some time making the wrapping from C++ into R (and back for that matter) a lot easier via our Rcpp package.

And because linear algebra is already such a well-understood and coded-for field, Armadillo, a current, modern, plesant, well-documted, small, templated, ... library was a very natural fit for our first extended wrapper: RcppArmadillo.

This has caught the attention of other MCMC users as well. I gave a one-day work at the U of Rochester business school last summer, and have help another researcher in the MidWest with similar explorations. Give RcppArmadillo a try -- it works well, is actively maintained (new Armadillo release 1.1.4 today, I will make a new RcppArmadillo later) and supported.

And because I just luuv this example so much, here is a quick "fast" version of lm() returning coefficient and std.errors:

extern "C" SEXP fastLm(SEXP ys, SEXP Xs) {

  try {
    Rcpp::NumericVector yr(ys);                 // creates Rcpp vector 
    Rcpp::NumericMatrix Xr(Xs);                 // creates Rcpp matrix 
    int n = Xr.nrow(), k = Xr.ncol();

    arma::mat X(Xr.begin(), n, k, false);       // avoids extra copy
    arma::colvec y(yr.begin(), yr.size(), false);

    arma::colvec coef = arma::solve(X, y);      // fit model y ~ X
    arma::colvec res = y - X*coef;              // residuals

    double s2 = std::inner_product(res.begin(), res.end(), 
                                   res.begin(), double())/(n - k);
                                            // std.errors of coefficients
    arma::colvec std_err = 
           arma::sqrt(s2 * arma::diagvec( arma::pinv(arma::trans(X)*X) ));  

    return Rcpp::List::create(Rcpp::Named("coefficients") = coef,
                              Rcpp::Named("stderr")       = std_err,
                              Rcpp::Named("df")           = n - k);

  } catch( std::exception &ex ) {
      forward_exception_to_r( ex );
  } catch(...) { 
      ::Rf_error( "c++ exception (unknown reason)" ); 
  }
  return R_NilValue; // -Wall
}

Lastly, you also get immediate prototyping via inline which may make 'time to code' faster.

Answer 2

I would strongly suggest that you have a look at RCpp and RcppArmadillo packages for R. Basically, you would not need to worry about the wrappers as they are already "included". Furthermore the syntactic sugar is really sweet (pun intended).

As a side remark, I would recommend that you have a look at JAGS, which does MCMC and its source code is in C++.

Answer 3

Boost Random from the Boost C++ libraries could be a good fit for you. In addition to many types of RNGs, it offers a variety of different distributions to draw from, such as

• Uniform (real)
• Uniform (unit sphere or arbitrary dimension)
• Bernoulli
• Binomial
• Cauchy
• Gamma
• Poisson
• Geometric
• Triangle
• Exponential
• Normal
• Lognormal

In addition, Boost Math complements the above distributions you can sample from with numerous density functions of many distributions. It also has several neat helper functions; just to give you an idea:

students_t dist(5);

cout << "CDF at t = 1 is " << cdf(dist, 1.0) << endl;
cout << "Complement of CDF at t = 1 is " << cdf(complement(dist, 1.0)) << endl;

for(double i = 10; i < 1e10; i *= 10)
{
   // Calculate the quantile for a 1 in i chance:
   double t = quantile(complement(dist, 1/i));
   // Print it out:
   cout << "Quantile of students-t with 5 degrees of freedom\n"
           "for a 1 in " << i << " chance is " << t << endl;
}

If you decided to use Boost, you also get to use its UBLAS library that features a variety of different matrix types and operations.

Answer 4

There are numerous C/C++ libraries out there, most focusing on a particular problem domain of (e.g. PDE solvers). There are two comprehensive libraries I can think of that you may find especially useful because they are written in C but have excellent Python wrappers already written.

1) IMSL C and PyIMSL

2) trilinos and pytrilinos

I have never used trilinos as the functionality is primarily on numerical analysis methods, but I use PyIMSL a lot for statistical work (and in a previous work life I developed the software too).

With respect to RNGs, here are the ones in C and Python in IMSL

DISCRETE

• random_binomial: Generates pseudorandom binomial numbers from a binomial distribution.
• random_geometric: Generates pseudorandom numbers from a geometric distribution.
• random_hypergeometric: Generates pseudorandom numbers from a hypergeometric distribution.
• random_logarithmic: Generates pseudorandom numbers from a logarithmic distribution.
• random_neg_binomial: Generates pseudorandom numbers from a negative binomial distribution.
• random_poisson: Generates pseudorandom numbers from a Poisson distribution.
• random_uniform_discrete: Generates pseudorandom numbers from a discrete uniform distribution.
• random_general_discrete: Generates pseudorandom numbers from a general discrete distribution using an alias method or optionally a table lookup method.

UNIVARIATE CONTINUOUS DISTRIBUTIONS

• random_beta: Generates pseudorandom numbers from a beta distribution.
• random_cauchy: Generates pseudorandom numbers from a Cauchy distribution.
• random_chi_squared: Generates pseudorandom numbers from a chi-squared distribution.
• random_exponential: Generates pseudorandom numbers from a standard exponential distribution.
• random_exponential_mix: Generates pseudorandom mixed numbers from a standard exponential distribution.
• random_gamma: Generates pseudorandom numbers from a standard gamma distribution.
• random_lognormal: Generates pseudorandom numbers from a lognormal distribution.
• random_normal: Generates pseudorandom numbers from a standard normal distribution.
• random_stable: Sets up a table to generate pseudorandom numbers from a general discrete distribution.
• random_student_t: Generates pseudorandom numbers from a Student's t distribution.
• random_triangular: Generates pseudorandom numbers from a triangular distribution.
• random_uniform: Generates pseudorandom numbers from a uniform (0, 1) distribution.
• random_von_mises: Generates pseudorandom numbers from a von Mises distribution.
• random_weibull: Generates pseudorandom numbers from a Weibull distribution.
• random_general_continuous: Generates pseudorandom numbers from a general continuous distribution.

MULTIVARIATE CONTINUOUS DISTRIBUTIONS

• random_normal_multivariate: Generates pseudorandom numbers from a multivariate normal distribution.
• random_orthogonal_matrix: Generates a pseudorandom orthogonal matrix or a correlation matrix.
• random_mvar_from_data: Generates pseudorandom numbers from a multivariate distribution determined from a given sample.
• random_multinomial: Generates pseudorandom numbers from a multinomial distribution.
• random_sphere: Generates pseudorandom points on a unit circle or K-dimensional sphere.
• random_table_twoway: Generates a pseudorandom two-way table.

ORDER STATISTICS

• random_order_normal: Generates pseudorandom order statistics from a standard normal distribution.
• random_order_uniform: Generates pseudorandom order statistics from a uniform (0, 1) distribution.

STOCHASTIC PROCESSES

• random_arma: Generates pseudorandom ARMA process numbers.
• random_npp: Generates pseudorandom numbers from a nonhomogeneous Poisson process.

SAMPLES AND PERMUTATIONS

• random_permutation: Generates a pseudorandom permutation.
• random_sample_indices: Generates a simple pseudorandom sample of indices.
• random_sample: Generates a simple pseudorandom sample from a finite population.

UTILITY FUNCTIONS

• random_option: Selects the uniform (0, 1) multiplicative congruential pseudorandom number generator.
• random_option_get: Retrieves the uniform (0, 1) multiplicative congruential pseudorandom number generator.
• random_seed_get: Retrieves the current value of the seed used in the IMSL random number generators.
• random_substream_seed_get: Retrieves a seed for the congruential generators that do not do shuffling that will generate random numbers beginning 100,000 numbers farther along.
• random_seed_set: Initializes a random seed for use in the IMSL random number generators.
• random_table_set: Sets the current table used in the shuffled generator.
• random_table_get: Retrieves the current table used in the shuffled generator.
• random_GFSR_table_set: Sets the current table used in the GFSR generator.
• random_GFSR_table_get: Retrieves the current table used in the GFSR generator.
• random_MT32_init: Initializes the 32-bit Mersenne Twister generator using an array.
• random_MT32_table_get: Retrieves the current table used in the 32-bit Mersenne Twister generator.
• random_MT32_table_set: Sets the current table used in the 32-bit Mersenne Twister generator.
• random_MT64_init: Initializes the 64-bit Mersenne Twister generator using an array.
• random_MT64_table_get: Retrieves the current table used in the 64-bit Mersenne Twister generator.
• random_MT64_table_set: Sets the current table used in the 64-bit Mersenne Twister generator.

LOW-DISCREPANCY SEQUENCE

• faure_next_point: Computes a shuffled Faure sequence.