Perpustakaan C ++ untuk komputasi statistik

Saya punya algoritma MCMC tertentu yang ingin saya porting ke C / C ++. Sebagian besar perhitungan mahal dalam C sudah melalui Cython, tapi saya ingin agar seluruh sampler ditulis dalam bahasa yang dikompilasi sehingga saya bisa menulis pembungkus untuk Python / R / Matlab / apa pun.

Setelah menyodok sekitar saya condong ke C ++. Beberapa perpustakaan yang relevan yang saya tahu adalah Armadillo (http://arma.sourceforge.net/) dan Scythe (http://scythe.wustl.edu/). Keduanya mencoba meniru beberapa aspek R / Matlab untuk mempermudah kurva belajar, yang sangat saya sukai. Sabit kotak sedikit lebih baik dengan apa yang ingin saya lakukan saya pikir. Secara khusus, RNG-nya mencakup banyak distribusi di mana Armadillo hanya memiliki seragam / normal, yang tidak nyaman. Armadillo tampaknya sedang dalam pengembangan yang cukup aktif sementara Scythe melihat rilis terakhirnya pada tahun 2007.

Jadi yang saya ingin tahu adalah apakah ada yang punya pengalaman dengan perpustakaan ini - atau orang lain yang hampir pasti saya lewatkan - dan jika demikian, apakah ada sesuatu untuk merekomendasikan satu di atas yang lain untuk ahli statistik yang sangat akrab dengan Python / R / Matlab tetapi kurang begitu dengan bahasa yang dikompilasi (tidak sepenuhnya bodoh, tetapi tidak benar-benar mahir ...).

Kami telah menghabiskan waktu membuat pembungkus dari C ++ menjadi R (dan kembali dalam hal ini) jauh lebih mudah melalui paket Rcpp kami .

Dan karena aljabar linier sudah menjadi bidang yang dipahami dengan baik dan dikodekan, Armadillo , arus, modern, plesant, didokumentasikan dengan baik, kecil, templated, ... perpustakaan sangat cocok untuk pembungkus pertama kami yang diperluas: RcppArmadillo .

Ini telah menarik perhatian pengguna MCMC lainnya juga. Saya memberikan satu hari kerja di sekolah bisnis U of Rochester musim panas lalu, dan telah membantu peneliti lain di MidWest dengan eksplorasi serupa. Berikan RcppArmadillo mencoba - bekerja dengan baik, secara aktif dipelihara (rilis Armadillo baru 1.1.4 hari ini, saya akan membuat RcppArmadillo baru nanti) dan didukung.

Dan karena saya sangat banyak memberikan contoh ini, berikut ini adalah versi cepat "cepat" untuk lm()koefisien kembali dan 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
}

Terakhir, Anda juga mendapatkan prototyping langsung melalui inline yang dapat membuat 'waktu untuk kode' lebih cepat.

Saya akan sangat menyarankan agar Anda melihat RCppdan RcppArmadillopaket untuk R. Pada dasarnya, Anda tidak perlu khawatir tentang pembungkus karena mereka sudah "disertakan". Selanjutnya gula sintaksis benar-benar manis (pun intended).

Sebagai komentar sampingan, saya akan merekomendasikan agar Anda melihat JAGS, yang tidak MCMC dan kode sumbernya di C ++.

Boost Random dari Boost C ++ library bisa sangat cocok untuk Anda. Selain banyak jenis RNG, ia menawarkan berbagai distribusi yang berbeda untuk diambil, seperti

• Seragam (nyata)
• Seragam (satuan bola atau dimensi arbitrer)
• Bernoulli
• Binomial
• Cauchy
• Gamma
• Poisson
• Geometris
• Segi tiga
• Eksponensial
• Normal
• Lognormal

Selain itu, Boost Math melengkapi distribusi di atas yang dapat Anda sampel dari dengan berbagai fungsi kepadatan dari banyak distribusi. Ini juga memiliki beberapa fungsi pembantu yang rapi; hanya untuk memberi Anda ide:

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

Jika Anda memutuskan untuk menggunakan Boost, Anda juga bisa menggunakan perpustakaan UBLAS-nya yang menampilkan berbagai jenis matriks dan operasi yang berbeda.

Ada banyak perpustakaan C / C ++ di luar sana, sebagian besar berfokus pada domain masalah tertentu (misalnya pemecah PDE). Ada dua perpustakaan komprehensif yang dapat saya pikirkan yang mungkin Anda temukan sangat berguna karena ditulis dalam C tetapi pembungkus Python yang sangat baik sudah ditulis.

1) IMSL C dan PyIMSL

2) trilinos dan pytrilinos

Saya tidak pernah menggunakan trilino karena fungsinya terutama pada metode analisis numerik, tetapi saya banyak menggunakan PyIMSL untuk pekerjaan statistik (dan dalam kehidupan kerja sebelumnya saya juga mengembangkan perangkat lunak).

Sehubungan dengan RNG, berikut adalah yang ada di C dan Python di IMSL

DISCRETE

• random_binomial: Menghasilkan angka binomial pseudorandom dari distribusi binomial.
• random_geometric: Menghasilkan angka pseudorandom dari distribusi geometris.
• random_hypergeometric: Menghasilkan angka pseudorandom dari distribusi hypergeometric.
• random_logarithmic: Menghasilkan angka pseudorandom dari distribusi logaritmik.
• random_neg_binomial: Menghasilkan angka pseudorandom dari distribusi binomial negatif.
• random_poisson: Menghasilkan angka pseudorandom dari distribusi Poisson.
• random_uniform_discrete: Menghasilkan angka pseudorandom dari distribusi seragam diskrit.
• random_general_discrete: Menghasilkan angka pseudorandom dari distribusi diskrit umum menggunakan metode alias atau opsional metode pencarian tabel.

DISTRIBUSI TERUS MENERUS

• random_beta: Menghasilkan angka pseudorandom dari distribusi beta.
• random_cauchy: Menghasilkan angka pseudorandom dari distribusi Cauchy.
• random_chi_squared: Menghasilkan angka pseudorandom dari distribusi chi-squared.
• random_exponential: Menghasilkan angka pseudorandom dari distribusi eksponensial standar.
• random_exponential_mix: Menghasilkan angka campuran pseudorandom dari distribusi eksponensial standar.
• random_gamma: Menghasilkan angka pseudorandom dari distribusi gamma standar.
• random_lognormal: Menghasilkan angka pseudorandom dari distribusi lognormal.
• random_normal: Menghasilkan angka pseudorandom dari distribusi normal standar.
• random_stable: Mengatur tabel untuk menghasilkan angka pseudorandom dari distribusi diskrit umum.
• random_student_t: Menghasilkan angka pseudorandom dari distribusi t Student.
• random_triangular: Menghasilkan angka pseudorandom dari distribusi segitiga.
• random_uniform: Menghasilkan angka pseudorandom dari distribusi seragam (0, 1).
• random_von_mises: Menghasilkan angka pseudorandom dari distribusi von Mises.
• random_weibull: Menghasilkan angka pseudorandom dari distribusi Weibull.
• random_general_continuous: Menghasilkan angka pseudorandom dari distribusi berkelanjutan umum.

DISTRIBUSI BERGANDA BERGANDA

• random_normal_multivariate: Menghasilkan angka pseudorandom dari distribusi normal multivarian.
• random_orthogonal_matrix: Menghasilkan matriks ortogonal pseudorandom atau matriks korelasi.
• random_mvar_from_data: Menghasilkan angka pseudorandom dari distribusi multivariat yang ditentukan dari sampel yang diberikan.
• random_multinomial: Menghasilkan angka pseudorandom dari distribusi multinomial.
• random_sphere: Menghasilkan titik pseudorandom pada lingkaran satuan atau bola dimensi K.
• random_table_twoway: Menghasilkan tabel dua arah pseudorandom.

STATISTIK PESANAN

• random_order_normal: Menghasilkan statistik pesanan pseudorandom dari distribusi normal standar.
• random_order_uniform: Menghasilkan statistik pesanan pseudorandom dari distribusi seragam (0, 1).

PROSES STOCHASTIC

• random_arma: Menghasilkan nomor proses ARMA pseudorandom.
• random_npp: Menghasilkan angka pseudorandom dari proses Poisson nonhomogen.

SAMPEL DAN PERMUTASI

• random_permutation: Menghasilkan permutasi pseudorandom.
• random_sample_indices: Menghasilkan sampel indeks pseudorandom sederhana.
• random_sample: Menghasilkan sampel pseudorandom sederhana dari populasi yang terbatas.

FUNGSI UTILITAS

• random_option: Memilih generator nomor pseudorandom multiplanatif seragam (0, 1).
• random_option_get: Mengambil generator nomor pseudorandom multipelatif seragam (0, 1).
• random_seed_get: Mengambil nilai saat ini dari benih yang digunakan dalam generator nomor acak IMSL.
• random_substream_seed_get: Mengambil seed untuk generator kongruensi yang tidak melakukan pengocokan yang akan menghasilkan bilangan acak yang dimulai dengan angka 100.000 lebih jauh.
• random_seed_set: Menginisialisasi benih acak untuk digunakan dalam generator nomor acak IMSL.
• random_table_set: Mengatur tabel saat ini digunakan dalam generator yang dikocok.
• random_table_get: Mengambil tabel saat ini yang digunakan dalam generator yang diacak.
• random_GFSR_table_set: Mengatur tabel saat ini yang digunakan dalam generator GFSR.
• random_GFSR_table_get: Mengambil tabel saat ini yang digunakan dalam generator GFSR.
• random_MT32_init: Menginisialisasi generator Mersenne Twister 32-bit menggunakan array.
• random_MT32_table_get: Mengambil tabel saat ini yang digunakan dalam generator Mersenne Twister 32-bit.
• random_MT32_table_set: Mengatur tabel saat ini yang digunakan dalam generator Mersenne Twister 32-bit.
• random_MT64_init: Menginisialisasi generator Mersenne Twister 64-bit menggunakan array.
• random_MT64_table_get: Mengambil tabel saat ini yang digunakan dalam generator Mersenne Twister 64-bit.
• random_MT64_table_set: Mengatur tabel saat ini yang digunakan dalam generator Mersenne Twister 64-bit.

URUTAN KECERDASAN RENDAH

• faure_next_point: Menghitung urutan Faure yang dikocok.