Bibliotecas C++ para el cálculo estadístico

Question

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 husmear me inclino por C++. Un par de bibliotecas relevantes que conozco son Armadillo (http://arma.sourceforge.net/) y Scythe (http://scythe.wustl.edu/). Ambas tratan de emular algunos aspectos de R/Matlab para facilitar la curva de aprendizaje, lo cual me gusta mucho. Scythe se ajusta un poco mejor a lo que quiero hacer, creo. En particular, su RNG incluye un montón de distribuciones donde Armadillo sólo tiene uniforme/normal, lo que es un inconveniente. Armadillo parece estar en desarrollo bastante activo mientras que Scythe vio su último lanzamiento en 2007.

Así que lo que me pregunto es si alguien tiene experiencia con estas librerías -o con otras que casi seguro he pasado por alto- y si es así, si hay algo que recomiende una sobre las otras para un estadístico muy familiarizado con Python/R/Matlab pero menos con los lenguajes compilados (no es completamente ignorante, pero no es exactamente competente...).

Answer 1

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

Y porque el álgebra lineal es ya un campo tan bien entendido y codificado, Armadillo Una librería actual, moderna, abundante, bien documentada, pequeña, con plantillas, ... era un ajuste muy natural para nuestra primera envoltura extendida: RcppArmadillo .

Esto ha llamado la atención de otros usuarios de MCMC. El verano pasado impartí una jornada de trabajo en la escuela de negocios de la Universidad de Rochester, y he ayudado a otro investigador del Medio Oeste con exploraciones similares. Regale RcppArmadillo una prueba -- funciona bien, se mantiene activamente (la nueva versión de Armadillo 1.1.4 de hoy, haré un nuevo RcppArmadillo más tarde) y tiene soporte.

Y como me encanta este ejemplo, aquí hay una versión rápida de lm() devolviendo el coeficiente y los 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 se obtiene la creación inmediata de prototipos a través de en línea lo que puede hacer que el "tiempo de codificación" sea más rápido.

Answer 2

Le sugiero encarecidamente que eche un vistazo a RCpp y RcppArmadillo paquetes para R . Básicamente, no tendría que preocuparse por las envolturas, ya que están "incluidas". Además, el azúcar sintáctico es realmente dulce (juego de palabras).

Como comentario adicional, le recomiendo que eche un vistazo a JAGS que hace MCMC y su código fuente está en C++.

Answer 3

Boost Random de las librerías Boost C++ podría ser una buena opción para ti. Además de muchos tipos de RNGs, ofrece una variedad de distribuciones diferentes de las que se puede sacar provecho, tales como

• Uniforme (real)
• Uniforme (esfera unitaria o dimensión arbitraria)
• Bernoulli
• Binomio
• Cauchy
• Gamma
• Poisson
• Geometría
• Triángulo
• Exponencial
• Normal
• Lognormal

Además, Impulsar las matemáticas complementa las distribuciones anteriores de las que puede tomar muestras con numerosas funciones de densidad de muchas distribuciones. También cuenta con varias funciones de ayuda muy útiles, para que te hagas 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 has decidido utilizar Boost, también podrás utilizar su biblioteca UBLAS que cuenta con una variedad de tipos de matrices y operaciones diferentes.

Answer 4

Existen numerosas bibliotecas de C/C++, la mayoría de ellas centradas en un dominio de problemas concreto (por ejemplo, solucionadores de EDP). Hay dos bibliotecas completas que se me ocurren que pueden ser 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 utilizado trilinos, ya que la funcionalidad es principalmente en los métodos de análisis numérico, pero yo uso PyIMSL mucho para el trabajo estadístico (y en una vida de trabajo anterior he desarrollado el software también).

Con respecto a los RNGs, aquí están 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 tablas.

DISTRIBUCIONES CONTINUAS UNIVARIANTES

• 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: Establece 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 MULTIVARIANTES

• random_normal_multivariate: Genera números pseudoaleatorios a partir de una distribución normal multivariante.
• matriz_ortogonal_aleatoria: 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 multivariable 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 en una esfera de K dimensiones.
• random_table_twoway: Genera una tabla pseudoaleatoria de dos vías.

ESTADÍSTICAS DE PEDIDOS

• 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 procesos ARMA pseudoaleatorios.
• random_npp: Genera números pseudoaleatorios a partir de un proceso de Poisson no homogéneo.

MUESTRAS Y PERMUTACIONES

• Permutación_aleatoria: Genera una permutación pseudoaleatoria.
• índices_de_muestra_aleatorios: Genera una muestra simple pseudoaleatoria de índices.
• muestra_aleatoria: Genera una muestra simple pseudoaleatoria de una población finita.

FUNCIONES DE UTILIDAD

• random_option: Selecciona el generador de números pseudoaleatorios congruente multiplicativo uniforme (0, 1).
• random_option_get: Recupera el generador de números pseudoaleatorios congruente 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 congruentes que no barajan que generarán números aleatorios empezando 100.000 números más adelante.
• random_seed_set: Inicializa una semilla aleatoria para su uso 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 de Mersenne Twister de 32 bits utilizando un array.
• random_MT32_table_get: Recupera la tabla actual utilizada en el generador de 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 de Mersenne Twister de 64 bits utilizando una matriz.
• random_MT64_table_get: Recupera la tabla actual utilizada en el generador de 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.