用於統計計算的 C++ 庫

我有一個特定的 MCMC 算法，我想將它移植到 C/C++。大部分昂貴的計算已經通過 Cython 用 C 語言完成，但我希望用編譯語言編寫整個採樣器，這樣我就可以為 Python/R/Matlab/whatever 編寫包裝器。

在閒逛之後，我傾向於 C++。我知道的幾個相關庫是 Armadillo (http://arma.sourceforge.net/) 和 Scythe (http://scythe.wustl.edu/)。兩者都試圖模仿 R/Matlab 的某些方面來簡化學習曲線，我非常喜歡這一點。我認為 Scythe 與我想做的事情相比要好一些。特別是，它的 RNG 包含很多分佈，其中犰狳只有統一/正態，這很不方便。Armadillo 似乎處於相當活躍的開發階段，而 Scythe 的最後一次發布是在 2007 年。

所以我想知道是否有人有使用這些庫的經驗——或者我幾乎肯定錯過的其他庫——如果是這樣，對於非常熟悉 Python/R/Matlab 的統計學家，是否有什麼可以推薦的但對於編譯語言則更少（並非完全無知，但並不完全精通……）。

通過我們的Rcpp包 ，我們花了一些時間使從 C++ 到R的包裝（以及返回）變得更加容易。

而且因為線性代數已經是一個很好理解和編碼的領域，犰狳，一個當前的、現代的、令人愉快的、記錄良好的、小型的、模板化的……庫非常適合我們的第一個擴展包裝器：RcppArmadillo .

這也引起了其他 MCMC 用戶的注意。去年夏天，我在羅切斯特大學商學院做了一天的工作，並幫助了中西部的另一位研究人員進行了類似的探索。試試 RcppArmadillo -它運行良好，得到積極維護（今天發布了新的 Armadillo 1.1.4，稍後我將製作一個新的 RcppArmadillo）並得到支持。

而且因為我只是對這個示例進行瞭如此多的操作，所以這裡有一個lm()返回係數和 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
}

最後，您還可以通過內聯獲得即時原型設計，這可能會使“編碼時間”更快。

引用自：https://stats.stackexchange.com/questions/7358