Biblioteki C ++ do obliczeń statystycznych

23

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 ...).

JMS
źródło

Odpowiedzi:

18

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”.

Dirk Eddelbuettel
źródło
Dzięki Dirk - miałem wrażenie, że odpowiesz wcześniej niż później :). Biorąc pod uwagę, że chcę kodu, który mogę wywoływać z innego oprogramowania (głównie Python, ale także Matlab), być może dobrym przepływem pracy byłoby prototypowanie w Rcpp / RcppArmadillo, a następnie przejście na „prosty” Armadillo? Składnia itp. Wygląda bardzo podobnie.
JMS
1
Mam nadzieję, że uznasz to za pomocne.
Dirk Eddelbuettel
Re twoje drugie pytanie z edycji: Jasne. Pancernik jest zależny od niewielkiego, lub w naszym przypadku, niczego poza R. Rcpp / RcppArmadillo nie pomógłby ci w interfejsie i testowaniu prototypowego kodu, który może być ponownie użyty samodzielnie lub z opakowaniami Python i Matlab, które możesz dodać później. Conrad może mieć wskaźniki do czegoś; Nie mam żadnych dla Pythona ani Matlaba.
Dirk Eddelbuettel
Przepraszam za wyciągnięcie dywanu :) Chcę, aby klawisz Enter zwrócił karetkę, ale zamiast tego przesyła mój komentarz. W każdym razie, dziękuję za pomoc - bawiłam się majsterkując i przeglądając dziś listę mailingową Rcpp.
JMS
8

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 ++.

zwiastun
źródło
2
Chciałbym to poprzeć. Jeśli szukasz szybki i łatwy sposób, aby interfejs skompilowany kod z R, Rcppze RcppArmadillojest droga. Edycja: Korzystając z Rcpp, masz również dostęp do wszystkich RNG zaimplementowanych w kodzie C leżącym u podstaw R.
fabians
Dziękuję za wotum zaufania. Już miałem to zaproponować ;-)
Dirk Eddelbuettel
7

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.

mavam
źródło
Dzięki za wskazówkę. Boost wygląda jak duży młotek dla mojego małego paznokcia, ale dojrzały i utrzymany.
JMS
Nie jestem pewien, czy boot :: math :: binomial_distribution ma taką samą funkcję, jak zaimplementowana w R binom.test () dwustronna. Zajrzałem do referencji i nie mogłem znaleźć tej funkcji. Próbowałem to zaimplementować i nie jest to trywialne!
Kemin Zhou
1

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.
Josh Hemann
źródło