Высокоточная сумма с R, включающая биномиальные коэффициенты и логарифмы

#r #precision #logarithm #binomial-coefficients

Вопрос:

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

Например:

 library(Rmpfr) binom lt;- function(n,i) {factorial(n)/(factorial(n-i)*factorial(i))}  i lt;- 30 n lt;-60 Ui lt;- rep(0,i) for (k in (0:(i-1))) {  Ui[k 1] lt;- (-1)^(i-1-k) * binom(i-1,k)/(n-k) * log(n-k) } U lt;- sum(mpfr(Ui, 1024))

возвращает 7.2395....e-10 , что очень далеко от фактического ответа, возвращенного Mathematica, то есть -5.11...e-20 .

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

Редактировать

Вот код Mathematica, который вычисляет ту же сумму. Он работает с целыми числами и преобразуется в реальные числа только после окончания суммы. Я увеличил количество сообщаемых десятичных знаков.

Причина этого?

В конце концов, мне нужно получить соотношение двух чисел, полученных с помощью аналогичных вычислений. Когда два числа отклоняются на пару порядков, полученное соотношение просто непредсказуемо.

1. начните с использования встроенного choose() выражения вместо вашего binom() выражения и посмотрите, поможет ли это.

2. @BenBolker: нет!

3. Да, я думаю, что это может быть непростой задачей. (1) Вы на 100% уверены, что нет ошибок/различий между тем, что вы реализовали здесь, и тем, что вы просили в Mathematica? (2) Я попробовал несколько вещей (упорядочение/группировка терминов по-разному при расчете), но пока безуспешно

4. Не могли бы вы, пожалуйста, показать свой код Mathematica?

5. @BenBolker Я добавил код Mathematica в правку.

Ответ №1:

Вам нужно работать с mpfr объектами в течение всего расчета, а не только при суммировании:

 library(Rmpfr)  i lt;- 30 n lt;- 60 k lt;- 0:(i - 1)  nk lt;- mpfr(n - k, 128) (U lt;- sum((-1)^(i-1-k)*choose(i-1,k)/(nk)*log(nk))) #gt; 1 'mpfr' number of precision 128 bits  #gt; [1] -5.110333215290518581300810256453669394729e-20  nk lt;- mpfr(n - k, 256) (U lt;- sum((-1)^(i-1-k)*choose(i-1,k)/(nk)*log(nk))) #gt; 1 'mpfr' number of precision 256 bits  #gt; [1] -5.110333215285320173235309727002720346864555872897902728222060861935229197560667e-20  nk lt;- mpfr(n - k, 512) (U lt;- sum((-1)^(i-1-k)*choose(i-1,k)/(nk)*log(nk))) #gt; 1 'mpfr' number of precision 512 bits  #gt; [1] -5.1103332152853201732353097270027203468645558728979134452318939958128833820370490135678222208577277855238767473116630391351888405531035522832949562601913591e-20