Bibliotecas C ++ para computación estadística

Tengo un algoritmo MCMC particular que me gustaría portar a C / C ++. Gran parte del cálculo costoso ya está en C a través de Cython, pero quiero tener todo el muestreador escrito en un lenguaje compilado para poder escribir envoltorios para Python / R / Matlab / lo que sea.

Después de hurgar, me estoy inclinando hacia C ++. Un par de bibliotecas relevantes que conozco son Armadillo (http://arma.sourceforge.net/) y Scythe (http://scythe.wustl.edu/). Ambos intentan emular algunos aspectos de R / Matlab para facilitar la curva de aprendizaje, lo que me gusta mucho. Scythe cuadra un poco mejor con lo que quiero hacer, creo. En particular, su RNG incluye muchas distribuciones donde Armadillo solo tiene uniforme / normal, lo cual es inconveniente. Armadillo parece estar bajo un desarrollo bastante activo, mientras que Scythe vio su último lanzamiento en 2007.

Entonces, lo que me pregunto es si alguien tiene experiencia con estas bibliotecas, u otras que seguramente he extrañado, y si es así, si hay algo para recomendar una sobre las otras para un estadístico muy familiarizado con Python / R / Matlab pero menos con lenguajes compilados (no completamente ignorantes, pero no exactamente competentes ...).

Hemos pasado algún tiempo haciendo que el ajuste de C ++ a R (y viceversa) sea mucho más fácil a través de nuestro paquete Rcpp .

Y debido a que el álgebra lineal ya es un campo tan bien entendido y codificado, Armadillo , una biblioteca actual, moderna, agradable, bien documentada, pequeña, con plantillas, ... era muy adecuada para nuestro primer envoltorio extendido: RcppArmadillo .

Esto también ha llamado la atención de otros usuarios de MCMC. Di un trabajo de un día en la escuela de negocios de la U de Rochester el verano pasado, y he ayudado a otro investigador en el Medio Oeste con exploraciones similares. Dar RcppArmadillo una oportunidad - que funciona bien, se mantiene activa (nueva versión 1.1.4 Armadillo hoy, voy a hacer una nueva RcppArmadillo más adelante) y apoyado.

Y debido a que me encanta este ejemplo, aquí hay una versión rápida "rápida" de lm()coeficientes de retorno y errores estándar:

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
}

Por último, también obtiene prototipos inmediatos a través de la línea, lo que puede hacer que el "tiempo de codificación" sea más rápido.

Sugeriría encarecidamente que eche un vistazo RCppy los RcppArmadillopaquetes para R. Básicamente, no necesitará preocuparse por los envoltorios ya que ya están "incluidos". Además, el azúcar sintáctico es realmente dulce (juego de palabras).

Como comentario adicional, recomendaría que eche un vistazo a JAGSMCMC y su código fuente está en C ++.

Boost Random de las bibliotecas Boost C ++ podría ser una buena opción para usted. Además de muchos tipos de RNG, ofrece una variedad de distribuciones diferentes para extraer, como

• Uniforme (real)
• Uniforme (esfera unitaria o dimensión arbitraria)
• Bernoulli
• Binomio
• Cauchy
• Gama
• Poisson
• Geométrico
• Triángulo
• Exponencial
• Normal
• Lognormal

Además, Boost Math complementa las distribuciones anteriores de las que puede muestrear con numerosas funciones de densidad de muchas distribuciones. También tiene varias funciones de ayuda ordenadas; solo para darte una 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;
}

Si decidió usar Boost, también puede usar su biblioteca UBLAS que presenta una variedad de diferentes tipos de matriz y operaciones.

Existen numerosas bibliotecas C / C ++, la mayoría centradas en un dominio de problema particular (por ejemplo, solucionadores de PDE). Se me ocurren dos bibliotecas completas que puede encontrar especialmente útiles porque están escritas en C pero tienen excelentes envoltorios de Python ya escritos.

1) IMSL C y PyIMSL

2) trilinos y pytrilinos

Nunca he usado trilinos ya que la funcionalidad se basa principalmente en métodos de análisis numérico, pero uso mucho PyIMSL para el trabajo estadístico (y en una vida laboral anterior también desarrollé el software).

Con respecto a los RNG, estos son los de C y Python en IMSL

DISCRETO

• random_binomial: genera números binomiales pseudoaleatorios a partir de una distribución binomial.
• random_geometric: genera números pseudoaleatorios a partir de una distribución geométrica.
• random_hypergeometric: genera números pseudoaleatorios a partir de una distribución hipergeométrica.
• random_logarithmic: genera números pseudoaleatorios a partir de una distribución logarítmica.
• random_neg_binomial: genera números pseudoaleatorios a partir de una distribución binomial negativa.
• random_poisson: genera números pseudoaleatorios a partir de una distribución de Poisson.
• random_uniform_discrete: genera números pseudoaleatorios a partir de una distribución uniforme discreta.
• random_general_discrete: genera números pseudoaleatorios a partir de una distribución discreta general utilizando un método de alias u opcionalmente un método de búsqueda de tabla.

DISTRIBUCIONES CONTINUAS UNIVARIADAS

• random_beta: genera números pseudoaleatorios a partir de una distribución beta.
• random_cauchy: genera números pseudoaleatorios a partir de una distribución Cauchy.
• random_chi_squared: genera números pseudoaleatorios a partir de una distribución chi-cuadrado.
• random_exponential: genera números pseudoaleatorios a partir de una distribución exponencial estándar.
• random_exponential_mix: genera números mixtos pseudoaleatorios a partir de una distribución exponencial estándar.
• random_gamma: genera números pseudoaleatorios a partir de una distribución gamma estándar.
• random_lognormal: genera números pseudoaleatorios a partir de una distribución lognormal.
• random_normal: genera números pseudoaleatorios a partir de una distribución normal estándar.
• random_stable: configura una tabla para generar números pseudoaleatorios a partir de una distribución discreta general.
• random_student_t: genera números pseudoaleatorios a partir de una distribución t de Student.
• random_triangular: genera números pseudoaleatorios a partir de una distribución triangular.
• random_uniform: genera números pseudoaleatorios a partir de una distribución uniforme (0, 1).
• random_von_mises: genera números pseudoaleatorios a partir de una distribución de von Mises.
• random_weibull: genera números pseudoaleatorios a partir de una distribución Weibull.
• random_general_continuous: genera números pseudoaleatorios a partir de una distribución continua general.

DISTRIBUCIONES CONTINUAS MULTIVARIAS

• random_normal_multivariate: genera números pseudoaleatorios a partir de una distribución normal multivariada.
• random_orthogonal_matrix: genera una matriz ortogonal pseudoaleatoria o una matriz de correlación.
• random_mvar_from_data: genera números pseudoaleatorios a partir de una distribución multivariada determinada a partir de una muestra dada.
• random_multinomial: genera números pseudoaleatorios a partir de una distribución multinomial.
• random_sphere: genera puntos pseudoaleatorios en un círculo unitario o esfera K-dimensional.
• random_table_twoway: genera una tabla bidireccional pseudoaleatoria.

ESTADÍSTICAS DE PEDIDO

• random_order_normal: genera estadísticas de orden pseudoaleatorio a partir de una distribución normal estándar.
• random_order_uniform: genera estadísticas de orden pseudoaleatorio a partir de una distribución uniforme (0, 1).

PROCESOS ESTOCÁSTICOS

• random_arma: genera números de proceso ARMA pseudoaleatorios.
• random_npp: genera números pseudoaleatorios a partir de un proceso de Poisson no homogéneo.

MUESTRAS Y PERMUTACIONES

• random_permutation: genera una permutación pseudoaleatoria.
• random_sample_indices: genera una muestra pseudoaleatoria simple de índices.
• random_sample: genera una muestra pseudoaleatoria simple a partir de una población finita.

FUNCIONES DE UTILIDAD

• random_option: selecciona el generador de números pseudoaleatorio congruencial multiplicativo uniforme (0, 1).
• random_option_get: recupera el generador de número pseudoaleatorio congruencial multiplicativo uniforme (0, 1).
• random_seed_get: recupera el valor actual de la semilla utilizada en los generadores de números aleatorios IMSL.
• random_substream_seed_get: recupera una semilla para los generadores congruenciales que no realizan barajado que generará números aleatorios que comenzarán 100,000 números más adelante.
• random_seed_set: Inicializa una semilla aleatoria para usar en los generadores de números aleatorios IMSL.
• random_table_set: establece la tabla actual utilizada en el generador aleatorio.
• random_table_get: recupera la tabla actual utilizada en el generador aleatorio.
• random_GFSR_table_set: establece la tabla actual utilizada en el generador GFSR.
• random_GFSR_table_get: recupera la tabla actual utilizada en el generador GFSR.
• random_MT32_init: Inicializa el generador Mersenne Twister de 32 bits utilizando una matriz.
• random_MT32_table_get: recupera la tabla actual utilizada en el generador Mersenne Twister de 32 bits.
• random_MT32_table_set: establece la tabla actual utilizada en el generador Mersenne Twister de 32 bits.
• random_MT64_init: Inicializa el generador Mersenne Twister de 64 bits utilizando una matriz.
• random_MT64_table_get: recupera la tabla actual utilizada en el generador Mersenne Twister de 64 bits.
• random_MT64_table_set: establece la tabla actual utilizada en el generador Mersenne Twister de 64 bits.

SECUENCIA DE BAJA DISCREPANCIA

• faure_next_point: calcula una secuencia de Faure barajada.