C ++ библиотеки для статистических вычислений

23

У меня есть определенный алгоритм MCMC, который я хотел бы перенести на C / C ++. Большая часть дорогостоящих вычислений выполняется на C уже через Cython, но я хочу, чтобы весь сэмплер был написан на скомпилированном языке, чтобы я мог просто писать оболочки для Python / R / Matlab / что угодно.

После осмотра я склоняюсь к C ++. Я знаю несколько соответствующих библиотек: Armadillo (http://arma.sourceforge.net/) и Scythe (http://scythe.wustl.edu/). Оба пытаются подражать некоторым аспектам R / Matlab, чтобы облегчить процесс обучения, что мне очень нравится. Я думаю, что косы немного лучше, чем я хочу. В частности, его RNG включает в себя множество дистрибутивов, где Armadillo имеет только равномерное / нормальное, что неудобно. Armadillo, кажется, находится в стадии довольно активной разработки, в то время как Scythe выпустила свой последний релиз в 2007 году.

Так что мне интересно, есть ли у кого-то опыт работы с этими библиотеками - или с другими, которые я почти наверняка пропустил - и если да, то есть ли что-то, что можно было бы порекомендовать одному из других для статистика, хорошо знакомого с Python / R / Matlab но не так с скомпилированными языками (не совсем невежественными, но не совсем опытными ...).

JMS
источник

Ответы:

18

Мы потратили некоторое время на то, чтобы упростить перенос из C ++ в R (и обратно) с помощью нашего пакета Rcpp .

И поскольку линейная алгебра уже является таким хорошо понятым и закодированным полем, Armadillo , современная, современная, приятная, хорошо документированная, небольшая, шаблонная, ... библиотека была очень естественной для нашей первой расширенной оболочки: RcppArmadillo ,

Это привлекло внимание и других пользователей MCMC. Прошлым летом я провел однодневную работу в бизнес-школе U of Rochester и помог другим исследователям в Среднем Западе провести аналогичные исследования. Дайте 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
}

И, наконец, вы также получаете немедленное прототипирование через inline, что может ускорить «кодирование».

Дирк Эддельбюттель
источник
Спасибо, Дирк - у меня было чувство, что ты ответишь скорее раньше, чем позже :). Учитывая, что мне нужен код, который я могу вызывать из другого программного обеспечения (в основном Python, но также и Matlab), возможно, хорошим рабочим процессом будет создание прототипа в Rcpp / RcppArmadillo, а затем переход к «прямому» Armadillo? Синтаксис и т. Д. Выглядит очень похоже.
JMS
1
Надеюсь, вы нашли это полезным.
Дирк Эддельбюттель
Ваш второй вопрос от редакции: Конечно. Armadillo мало зависит, или, в нашем случае, ничто, кроме R. Rcpp / RcppArmadillo, не поможет вам в интерфейсе и тестировании прототипированного кода, который можно использовать повторно отдельно или с обертками Python и Matlab, которые вы можете добавить позже. Конрад может иметь указатели на что-то; У меня нет ни одного для Python или Matlab.
Дирк Эддельбюттель
Извините, что вытащил коврик :) Я хочу, чтобы клавиша ввода выдала возврат каретки, но вместо этого он отправляет мой комментарий. В любом случае, спасибо за вашу помощь - я наслаждался тем, что порылся и копался в списке рассылки Rcpp сегодня весь день.
JMS
8

Я бы настоятельно рекомендую вам взглянуть на RCppи RcppArmadilloпакеты для R. По сути, вам не нужно беспокоиться об оболочках, так как они уже включены. Кроме того, синтаксический сахар действительно сладкий (каламбур).

В качестве побочного замечания я бы порекомендовал вам взглянуть JAGS, что делает MCMC и его исходный код находится на C ++.

Тевкр
источник
2
Я хотел бы поддержать это. Если вы ищете быстрый и простой способ интерфейс скомпилированный код с R, Rcppс RcppArmadilloэто путь. Редактировать: Используя Rcpp, у вас также есть доступ ко всем ГСЧ, включенным в C-код, лежащий в основе R.
Фабианс
Спасибо за вотум доверия. Я собирался предложить то же самое ;-)
Dirk Eddelbuettel
7

Boost Random из библиотек Boost C ++ может подойти вам. В дополнение ко многим типам RNG, он предлагает множество различных распределений, таких как:

  • Униформа (реальная)
  • Униформа (единичная сфера или произвольное измерение)
  • Бернулли
  • бином
  • коши
  • Гамма
  • Пуассон
  • геометрический
  • Треугольник
  • экспоненциальный
  • Обычный
  • Логарифмическая

Кроме того, Boost Math дополняет вышеприведенные распределения, из которых можно выполнить выборку, множеством функций плотности многих распределений. Он также имеет несколько полезных вспомогательных функций; просто чтобы дать вам представление:

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

Если вы решили использовать Boost, вы также можете использовать его библиотеку UBLAS, которая содержит множество различных типов матриц и операций.

mavam
источник
Спасибо за совет. Boost выглядит как большой молоток для моего маленького ногтя, но зрелый и ухоженный.
JMS
Я не уверен, что boot :: math :: binomial_distribution имеет ту же функцию, что и реализованная в R binom.test (), двусторонняя. Я посмотрел в ссылку и не смог найти эту функцию. Я пытался реализовать это, и это не тривиально!
Кемин Чжоу
1

Существует множество библиотек C / C ++, большинство из которых сосредоточено на конкретной проблемной области (например, решатели PDE). Я могу подумать, что есть две всеобъемлющие библиотеки, которые могут оказаться особенно полезными, поскольку они написаны на C, но уже имеют отличные оболочки Python.

1) IMSL C и PyIMSL

2) трилино и питрилино

Я никогда не использовал трилино, так как функциональность в основном связана с методами численного анализа, но я часто использую PyIMSL для статистической работы (и в предыдущей жизни я тоже разрабатывал программное обеспечение).

Что касается RNG, вот те, в C и Python в IMSL

DISCRETE

  • random_binomial: генерирует псевдослучайные биномиальные числа из биномиального распределения.
  • random_geometric: генерирует псевдослучайные числа из геометрического распределения.
  • random_hypergeometric: генерирует псевдослучайные числа из гипергеометрического распределения.
  • random_logarithmic: генерирует псевдослучайные числа из логарифмического распределения.
  • random_neg_binomial: генерирует псевдослучайные числа из отрицательного биномиального распределения.
  • random_poisson: генерирует псевдослучайные числа из распределения Пуассона.
  • random_uniform_discrete: генерирует псевдослучайные числа из дискретного равномерного распределения.
  • random_general_discrete: генерирует псевдослучайные числа из общего дискретного распределения, используя метод псевдонимов или, возможно, метод поиска в таблице.

УНИВЕРСАЛЬНЫЕ НЕПРЕРЫВНЫЕ РАСПРЕДЕЛЕНИЯ

  • random_beta: генерирует псевдослучайные числа из бета-распределения.
  • random_cauchy: генерирует псевдослучайные числа из распределения Коши.
  • random_chi_squared: генерирует псевдослучайные числа из распределения хи-квадрат.
  • random_exponential: генерирует псевдослучайные числа из стандартного экспоненциального распределения.
  • random_exponential_mix: генерирует псевдослучайные смешанные числа из стандартного экспоненциального распределения.
  • random_gamma: генерирует псевдослучайные числа из стандартного гамма-распределения.
  • random_lognormal: генерирует псевдослучайные числа из логнормального распределения.
  • random_normal: генерирует псевдослучайные числа из стандартного нормального распределения.
  • random_stable: устанавливает таблицу для генерации псевдослучайных чисел из общего дискретного распределения.
  • random_student_t: генерирует псевдослучайные числа из t-распределения Стьюдента.
  • random_triangular: генерирует псевдослучайные числа из треугольного распределения.
  • random_uniform: генерирует псевдослучайные числа из равномерного (0, 1) распределения.
  • random_von_mises: генерирует псевдослучайные числа из распределения фон Мизеса.
  • random_weibull: генерирует псевдослучайные числа из распределения Вейбулла.
  • random_general_continuous: генерирует псевдослучайные числа из общего непрерывного распределения.

МНОГОМЕРНЫЕ НЕПРЕРЫВНЫЕ РАСПРЕДЕЛЕНИЯ

  • random_normal_multivariate: генерирует псевдослучайные числа из многомерного нормального распределения.
  • random_orthogonal_matrix: генерирует псевдослучайную ортогональную матрицу или матрицу корреляции.
  • random_mvar_from_data: генерирует псевдослучайные числа из многомерного распределения, определенного из заданной выборки.
  • random_multinomial: генерирует псевдослучайные числа из полиномиального распределения.
  • random_sphere: генерирует псевдослучайные точки на единичной окружности или K-мерной сфере.
  • random_table_twoway: генерирует псевдослучайную двустороннюю таблицу.

СТАТИСТИКА ЗАКАЗА

  • random_order_normal: генерирует статистику псевдослучайного порядка из стандартного нормального распределения.
  • random_order_uniform: генерирует статистику псевдослучайного порядка из равномерного (0, 1) распределения.

СТОХАСТИЧЕСКИЕ ПРОЦЕССЫ

  • random_arma: генерирует псевдослучайные номера процессов ARMA.
  • random_npp: генерирует псевдослучайные числа из неоднородного пуассоновского процесса.

ОБРАЗЦЫ И ГЕРМЕТИКИ

  • random_permutation: генерирует псевдослучайную перестановку.
  • random_sample_indices: генерирует простую псевдослучайную выборку индексов.
  • random_sample: генерирует простую псевдослучайную выборку из конечной совокупности.

ПОЛЕЗНЫЕ ФУНКЦИИ

  • random_option: выбирает равномерный (0, 1) мультипликативный конгруэнтный генератор псевдослучайных чисел.
  • random_option_get: извлекает равномерный (0, 1) мультипликативный конгруэнтный генератор псевдослучайных чисел.
  • random_seed_get: извлекает текущее значение начального числа, используемого в генераторах случайных чисел IMSL.
  • random_substream_seed_get: извлекает начальное число для конгруэнтных генераторов, которые не выполняют перемешивание, которые будут генерировать случайные числа, начинающиеся на 100 000 чисел дальше.
  • random_seed_set: Инициализирует случайное начальное число для использования в генераторах случайных чисел IMSL.
  • random_table_set: устанавливает текущую таблицу, используемую в перемешанном генераторе.
  • random_table_get: извлекает текущую таблицу, используемую в перемешанном генераторе.
  • random_GFSR_table_set: устанавливает текущую таблицу, используемую в генераторе GFSR.
  • random_GFSR_table_get: извлекает текущую таблицу, используемую в генераторе GFSR.
  • random_MT32_init: Инициализирует 32-битный генератор Mersenne Twister, используя массив.
  • random_MT32_table_get: извлекает текущую таблицу, используемую в 32-битном генераторе Mersenne Twister.
  • random_MT32_table_set: Устанавливает текущую таблицу, используемую в 32-битном генераторе Mersenne Twister.
  • random_MT64_init: Инициализирует 64-битный генератор Mersenne Twister, используя массив.
  • random_MT64_table_get: извлекает текущую таблицу, используемую в 64-битном генераторе Mersenne Twister.
  • random_MT64_table_set: устанавливает текущую таблицу, используемую в 64-битном генераторе Mersenne Twister.

НИЗКАЯ РАСПРОСТРАНЕННОСТЬ

  • faure_next_point: вычисляет перемешанную последовательность Faure.
Джош Хеманн
источник