Biblioteki C ++ do obliczeń statystycznych

Mam określony algorytm MCMC, który chciałbym przenieść do C / C ++. Wiele kosztownych obliczeń jest już napisanych w C przez Cython, ale chcę mieć cały sampler napisany w skompilowanym języku, aby móc po prostu pisać opakowania dla Python / R / Matlab / cokolwiek.

Po przeszukiwaniu skłaniam się ku C ++. Kilka odpowiednich bibliotek, które znam, to Armadillo (http://arma.sourceforge.net/) i Scythe (http://scythe.wustl.edu/). Oba starają się naśladować niektóre aspekty R / Matlab, aby ułatwić krzywą uczenia się, co bardzo mi się podoba. Myślę, że kosa trochę lepiej łączy się z tym, co chcę zrobić. W szczególności jego RNG obejmuje wiele dystrybucji, w których Armadillo ma tylko jednolity / normalny, co jest niewygodne. Wydaje się, że Armadillo jest dość aktywnie rozwijany, podczas gdy Scythe widział swoje ostatnie wydanie w 2007.

Zastanawiam się więc, czy ktoś ma doświadczenie z tymi bibliotekami - lub z innymi, które prawie na pewno tęskniłem - a jeśli tak, to czy jest coś, co poleciłoby je innym dla statystyki dobrze znającej Python / R / Matlab ale w mniejszym stopniu ze skompilowanymi językami (nie do końca ignoranccy, ale nie do końca biegli ...).

Poświęciliśmy trochę czasu na znacznie łatwiejsze pakowanie z C ++ do R (iz powrotem) dzięki naszemu pakietowi Rcpp .

A ponieważ algebra liniowa jest już tak dobrze zrozumiałym i zakodowanym polem, Armadillo , aktualna, nowoczesna, przyjemna, dobrze udokumentowana, mała, wzorowana ... biblioteka była bardzo naturalnym dopasowaniem do naszego pierwszego rozszerzonego opakowania: RcppArmadillo .

Zwróciło to również uwagę innych użytkowników MCMC. Zeszłego lata dałem jednodniową pracę w szkole biznesu U w Rochester i pomogłem innemu badaczowi z MidWest przy podobnych poszukiwaniach. Wypróbuj RcppArmadillo - działa dobrze, jest aktywnie utrzymywany (nowa wersja Armadillo 1.1.4 dzisiaj, zrobię nowy RcppArmadillo później) i obsługiwany.

A ponieważ tak bardzo opisałem ten przykład, oto szybka „szybka” wersja lm()zwracanego współczynnika i błędów standardowych:

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
}

Wreszcie otrzymujesz również natychmiastowe prototypowanie za pomocą wbudowanego narzędzia, co może przyspieszyć „czas kodowania”.

Zdecydowanie sugeruję, abyś spojrzał na RCppi RcppArmadillopakiety R. Zasadniczo nie trzeba się martwić o opakowania, ponieważ są one już „dołączone”. Ponadto cukier syntaktyczny jest naprawdę słodki (zamierzone słowo).

Na marginesie, polecam przyjrzeć się JAGS, co robi MCMC i jego kod źródłowy jest w C ++.

Boost Random z bibliotek Boost C ++ może być dla Ciebie dobrym rozwiązaniem. Oprócz wielu rodzajów RNG oferuje wiele różnych dystrybucji, z których można czerpać, takich jak

• Jednolity (prawdziwy)
• Jednolity (sfera jednostkowa lub dowolny wymiar)
• Bernoulli
• Dwumianowy
• Cauchy
• Gamma
• Poisson
• Geometryczny
• Trójkąt
• Wykładniczy
• Normalna
• Lognormal

Ponadto, Boost Math uzupełnia powyższe rozkłady, z których można próbkować, dzięki licznym funkcjom gęstości wielu rozkładów. Ma także kilka schludnych funkcji pomocniczych; aby dać ci pomysł:

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;
}

Jeśli zdecydujesz się użyć Boost, możesz także skorzystać z biblioteki UBLAS, która zawiera wiele różnych typów macierzy i operacji.

Istnieje wiele bibliotek C / C ++, z których większość koncentruje się na konkretnej dziedzinie problemów (np. Rozwiązaniach PDE). Są dwie obszerne biblioteki, o których myślę, że mogą być szczególnie przydatne, ponieważ są napisane w języku C, ale mają już doskonałe opakowanie w języku Python.

1) IMSL C i PyIMSL

2) trilinos i pytrilinos

Nigdy nie korzystałem z trilinos, ponieważ funkcjonalność dotyczy przede wszystkim metod analizy numerycznej, ale PyIMSL bardzo często używam do pracy statystycznej (w poprzednim życiu też opracowałem oprogramowanie).

W odniesieniu do RNG, oto te w C i Python w IMSL

ODDZIELNY

• random_binomial: Generuje pseudolosowe liczby dwumianowe z rozkładu dwumianowego.
• random_geometric: Generuje liczby pseudolosowe z rozkładu geometrycznego.
• random_hypergeometric: Generuje liczby pseudolosowe z rozkładu hipergeometrycznego.
• random_logarithmic: Generuje liczby pseudolosowe z rozkładu logarytmicznego.
• random_neg_binomial: Generuje liczby pseudolosowe z ujemnego rozkładu dwumianowego.
• random_poisson: Generuje liczby pseudolosowe z rozkładu Poissona.
• random_uniform_discrete: Generuje liczby pseudolosowe z dyskretnego rozkładu równomiernego.
• random_general_discrete: Generuje liczby pseudolosowe z ogólnej dystrybucji dyskretnej przy użyciu metody aliasowej lub opcjonalnie metody wyszukiwania w tabeli.

UNIWERSALNE CIĄGŁE DYSTRYBUCJE

• random_beta: Generuje pseudolosowe liczby z rozkładu beta.
• random_cauchy: Generuje liczby pseudolosowe z rozkładu Cauchy'ego.
• random_chi_squared: Generuje liczby pseudolosowe z rozkładu chi-kwadrat.
• random_exponential: Generuje liczby pseudolosowe ze standardowego rozkładu wykładniczego.
• random_exponential_mix: Generuje pseudolosowe liczby mieszane ze standardowego rozkładu wykładniczego.
• random_gamma: Generuje liczby pseudolosowe ze standardowego rozkładu gamma.
• random_lognormal: Generuje liczby pseudolosowe z rozkładu logarytmicznego.
• random_normal: Generuje liczby pseudolosowe ze standardowego rozkładu normalnego.
• random_stable: tworzy tabelę do generowania liczb pseudolosowych z ogólnej dystrybucji dyskretnej.
• random_student_t: Generuje liczby pseudolosowe z rozkładu t Studenta.
• random_triangular: Generuje liczby pseudolosowe z rozkładu trójkątnego.
• random_uniform: Generuje liczby pseudolosowe z jednolitego rozkładu (0, 1).
• random_von_mises: Generuje liczby pseudolosowe z rozkładu von Misesa.
• random_weibull: Generuje liczby pseudolosowe z rozkładu Weibulla.
• random_general_continuous: Generuje liczby pseudolosowe z ogólnego ciągłego rozkładu.

WIELOKROTNE CIĄGŁE DYSTRYBUCJE

• random_normal_multivariate: Generuje liczby pseudolosowe z wielowymiarowego rozkładu normalnego.
• random_orthogonal_matrix: Generuje pseudolosową macierz ortogonalną lub macierz korelacji.
• random_mvar_from_data: Generuje liczby pseudolosowe z rozkładu wielowymiarowego określonego z danej próbki.
• random_multinomial: Generuje liczby pseudolosowe z rozkładu wielomianowego.
• random_sphere: Generuje pseudolosowe punkty na okręgu jednostkowym lub kuli K-wymiarowej.
• random_table_twoway: Generuje pseudolosową dwukierunkową tabelę.

STATYSTYKA ZAMÓWIEŃ

• random_order_normal: Generuje statystyki kolejności pseudolosowej ze standardowego rozkładu normalnego.
• random_order_uniform: Generuje statystyki rzędu pseudolosowego z jednolitego rozkładu (0, 1).

PROCESY STOCHASTYCZNE

• random_arma: Generuje pseudolosowe numery procesów ARMA.
• random_npp: Generuje liczby pseudolosowe z niejednorodnego procesu Poissona.

PRÓBKI I UPRAWNIENIA

• random_permutation: Generuje pseudolosową permutację.
• random_sample_indices: Generuje prostą pseudolosową próbkę indeksów.
• random_sample: generuje prostą próbkę pseudolosową ze skończonej populacji.

FUNKCJE UŻYTKOWE

• opcja_ losowa: Wybiera jednolity (0, 1) multiplikatywny generator pseudolosowych liczb zbieżnych.
• random_option_get: Pobiera jednolity (0, 1) multiplikatywny generator pseudolosowych liczb zbieżnych.
• random_seed_get: pobiera bieżącą wartość zarodka używanego w generatorach liczb losowych IMSL.
• random_substream_seed_get: Pobiera ziarno dla generatorów kongruencjalnych, które nie wykonują tasowania, generując liczby losowe rozpoczynające się o 100 000 liczb dalej.
• random_seed_set: Inicjuje losowe ziarno do użycia w generatorach liczb losowych IMSL.
• random_table_set: Ustawia bieżącą tabelę używaną w losowym generatorze.
• random_table_get: pobiera bieżącą tabelę używaną w losowym generatorze.
• random_GFSR_table_set: Ustawia bieżącą tabelę używaną w generatorze GFSR.
• random_GFSR_table_get: Pobiera bieżącą tabelę używaną w generatorze GFSR.
• random_MT32_init: Inicjuje 32-bitowy generator Mersenne Twister przy użyciu tablicy.
• random_MT32_table_get: pobiera bieżącą tabelę używaną w 32-bitowym generatorze Mersenne Twister.
• random_MT32_table_set: Ustawia bieżącą tabelę używaną w 32-bitowym generatorze Mersenne Twister.
• random_MT64_init: Inicjuje 64-bitowy generator Mersenne Twister przy użyciu tablicy.
• random_MT64_table_get: pobiera bieżącą tabelę używaną w 64-bitowym generatorze Mersenne Twister.
• random_MT64_table_set: Ustawia bieżącą tabelę używaną w 64-bitowym generatorze Mersenne Twister.

SEKWENCJA O NISKIEJ ROZKŁADU

• faure_next_point: Oblicza potasowaną sekwencję Faure.