#include <boost/math/special_functions/beta.hpp>

namespace boost{ namespace math{
template <class T1, class T2, class T3>
calculated-result-type ibeta_inv(T1 a, T2 b, T3 p);
template <class T1, class T2, class T3, class Policy>
calculated-result-type ibeta_inv(T1 a, T2 b, T3 p, const Policy&);
template <class T1, class T2, class T3, class T4>
calculated-result-type ibeta_inv(T1 a, T2 b, T3 p, T4* py);
template <class T1, class T2, class T3, class T4, class Policy>
calculated-result-type ibeta_inv(T1 a, T2 b, T3 p, T4* py, const Policy&);
template <class T1, class T2, class T3>
calculated-result-type ibetac_inv(T1 a, T2 b, T3 q);
template <class T1, class T2, class T3, class Policy>
calculated-result-type ibetac_inv(T1 a, T2 b, T3 q, const Policy&);
template <class T1, class T2, class T3, class T4>
calculated-result-type ibetac_inv(T1 a, T2 b, T3 q, T4* py);
template <class T1, class T2, class T3, class T4, class Policy>
calculated-result-type ibetac_inv(T1 a, T2 b, T3 q, T4* py, const Policy&);
template <class T1, class T2, class T3>
calculated-result-type ibeta_inva(T1 b, T2 x, T3 p);
template <class T1, class T2, class T3, class Policy>
calculated-result-type ibeta_inva(T1 b, T2 x, T3 p, const Policy&);
template <class T1, class T2, class T3>
calculated-result-type ibetac_inva(T1 b, T2 x, T3 q);
template <class T1, class T2, class T3, class Policy>
calculated-result-type ibetac_inva(T1 b, T2 x, T3 q, const Policy&);
template <class T1, class T2, class T3>
calculated-result-type ibeta_invb(T1 a, T2 x, T3 p);
template <class T1, class T2, class T3, class Policy>
calculated-result-type ibeta_invb(T1 a, T2 x, T3 p, const Policy&);
template <class T1, class T2, class T3>
calculated-result-type ibetac_invb(T1 a, T2 x, T3 q);
template <class T1, class T2, class T3, class Policy>
calculated-result-type ibetac_invb(T1 a, T2 x, T3 q, const Policy&);
}} // namespaces

Description

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

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

Tip

Когда люди обычно говорят об обратной неполной бета-функции, они говорят об инвертировании по параметруx. Они реализованы здесь как ibeta_inv и ibetac_inv, и на сегодняшний день являются наиболее эффективными из обратных, представленных здесь. Обратные параметры параметровaиbнаходят применение в некоторых статистических приложениях, но должны быть вычислены с помощью довольно грубых численных методов и, следовательно, в несколько раз медленнее. Они реализованы здесь как ibeta_inva и ibeta_invb, а также дополняют версии ibetac_inva и ibetac_invb.

Тип возврата этих функций вычисляется с помощью.правила расчета типа результатапри вызове аргументов Т1... TN различных типов.

template <class T1, class T2, class T3>
calculated-result-type ibeta_inv(T1 a, T2 b, T3 p);
template <class T1, class T2, class T3, class Policy>
calculated-result-type ibeta_inv(T1 a, T2 b, T3 p, const Policy&);
template <class T1, class T2, class T3, class T4>
calculated-result-type ibeta_inv(T1 a, T2 b, T3 p, T4* py);
template <class T1, class T2, class T3, class T4, class Policy>
calculated-result-type ibeta_inv(T1 a, T2 b, T3 p, T4* py, const Policy&);

Возвращает значениеxтак, что:<p =ibeta(a, b,x);>и устанавливает<*py =1-x>, когда<py>параметр обеспечен и не является нулевым. Обратите внимание, что внутренне эта функция вычисляет в зависимости от того, что меньше<x>и<1-x>, и, следовательно, значение, присвоенное<*py>, не содержит ошибок отмены. Это означает, что даже если функция возвращает<1>, значение, сохраненное в<*py>, может быть ненулевым, хотя и очень маленьким.

Требуется:a,b >0и0<= p<= 1.

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

template <class T1, class T2, class T3>
calculated-result-type ibetac_inv(T1 a, T2 b, T3 q);
template <class T1, class T2, class T3, class Policy>
calculated-result-type ibetac_inv(T1 a, T2 b, T3 q, const Policy&);
template <class T1, class T2, class T3, class T4>
calculated-result-type ibetac_inv(T1 a, T2 b, T3 q, T4* py);
template <class T1, class T2, class T3, class T4, class Policy>
calculated-result-type ibetac_inv(T1 a, T2 b, T3 q, T4* py, const Policy&);

Возвращает значениеxтак, что:<q =ibetac(a, b,x);>и устанавливает<*py =1-x>, когда параметр<py>предусмотрен и не является нулевым. Обратите внимание, что внутренне эта функция вычисляет в зависимости от того, что меньше<x>и<1-x>, и, следовательно, значение, присвоенное<*py>, не содержит ошибок отмены. Это означает, что даже если функция возвращается<1>, значение, сохраненное в<*py>, может быть ненулевым, хотя и очень маленьким.

Требуется:a,b >0и0<= q<= 1.

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

template <class T1, class T2, class T3>
calculated-result-type ibeta_inva(T1 b, T2 x, T3 p);
template <class T1, class T2, class T3, class Policy>
calculated-result-type ibeta_inva(T1 b, T2 x, T3 p, const Policy&);

Возвращает значениетак, что:<p =ibeta(a, b,x);>

Требуется:b >0,0< x< 1и0<= p<= 1.

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

template <class T1, class T2, class T3>
calculated-result-type ibetac_inva(T1 b, T2 x, T3 p);
template <class T1, class T2, class T3, class Policy>
calculated-result-type ibetac_inva(T1 b, T2 x, T3 p, const Policy&);

Возвращает значениеaтак, что:<q =ibetac(a, b,x);>

Требуется:b >0,0< x< 1и0<= q<= 1.

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

template <class T1, class T2, class T3>
calculated-result-type ibeta_invb(T1 b, T2 x, T3 p);
template <class T1, class T2, class T3, class Policy>
calculated-result-type ibeta_invb(T1 b, T2 x, T3 p, const Policy&);

Возвращает значениеbтак, что:<p =ibeta(a, b,x);>

Требуется:a >0,0< x< 1и0<= p<= 1.

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

template <class T1, class T2, class T3>
calculated-result-type ibetac_invb(T1 b, T2 x, T3 p);
template <class T1, class T2, class T3, class Policy>
calculated-result-type ibetac_invb(T1 b, T2 x, T3 p, const Policy&);

Возвращает значениеbтак, что:<q =ibetac(a, b,x);>

Требуется:a >0,0< x< 1и0<= q<= 1.

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

Accuracy

Точность этих функций должна точно соответствовать точности регулярных неполных бета-функций. Однако обратите внимание, что в некоторых частях их домена эти функции могут быть чрезвычайно чувствительны к изменениям ввода, особенно когда аргументp(или его дополнениеq) очень близок к<0>или<1>.

Сравнения с другими библиотеками приведены ниже, обратите внимание, что наши тестовые данные используют некоторые довольно крайние случаи в неполной бета-функции, с которой многие другие библиотеки не справляются:

Testing

Существует два набора тестов:

Implementation of ibeta_inv and ibetac_inv

Эти две функции имеют общую реализацию.

Сначала вычисляется начальное приближение к x, затем последние несколько бит очищаются с помощью. Итерация Галлея. Предел итерации устанавливается до 1/2 числа битов в T, что является достаточным для того, чтобы гарантировать, что инверсии по меньшей мере столь же точны, как и обычные неполные бета-функции. В крайних случаях может потребоваться до 5 итераций, хотя обычно требуется только одна или две. Кроме того, количество требуемых итераций уменьшается с увеличениемaиb(которые обычно образуют более важные варианты использования).

Первоначальные предположения, используемые для итерации, получаются следующим образом:

Во-первых, напомним, что:

Мы можем начать с p или q и вычислить x или y. Кроме того, на любом этапе мы можем обменять a на b, p на q и x на y, если это приведет к более управляемой проблеме.

Для<a+b>=5>исходное предположение вычисляется с использованием способов, описанных в:

Асимптотическая инверсия неполной бета-функции, Н. М.Темме. Журнал вычислительной и прикладной математики 41 (1992) 145-157.

Почти симметричный корпус (раздел 2) используется для

В первую очередь необходимо решить функцию обратной ошибки. Метод точен по меньшей мере до 2 десятичных цифр, когда<a = 5>поднимается до по меньшей мере 8 цифр, когда<a = 105>.

Случай общей функции ошибки (раздел 3 статьи) используется для

и снова выражает обратную неполную бета в терминах обратной функции ошибки. Метод точен по меньшей мере до 2 десятичных цифр, когда<a+b = 5>поднимается до 11 цифр, когда<a+b = 105>. Однако, когда результат, как ожидается, будет очень маленьким, и когда a+b также мал, то его точность отпадает, в этом случае, когда p^1/a< 0,0025 тогда лучше использовать следующее в качестве первоначальной оценки:

Наконец, для всех других случаев, где<a+b> 5>используется метод раздела 4 статьи. Это выражает обратную неполную бета с точки зрения обратной неполной гамма-функции и, следовательно, значительно дороже, чем в других случаях. Однако метод точен по меньшей мере до 3 десятичных цифр, когда<a = 5>поднимается до по меньшей мере 10 цифр, когда<a = 105>. Этот метод ограничен a >b, и поэтому нам нужно произвести обмен a на b, p на q и x на y, когда это не так. Кроме того, когда p близко к 1, метод является неточным, если мы действительно хотим y, а не x в качестве вывода. Поэтому, когда q мало (<q1/p< 10-3>), мы используем:

Это дешевле, чем полный метод, и более точная оценка на q.

Когда a и b малы, в литературе явно не хватает информации о том, как действовать. Я очень благодарен профессору Нико Темме, который предоставил следующую информацию с большим терпением и объяснением с его стороны. Все последующие ошибки полностью принадлежат мне, а не профессору Темме.

Если a и b меньше 1, то в точке неполной бета-версии имеется точка перегиба<xs =(1-a)/(2-a -b)>. Поэтому, если<p > Ix(a,b)>мы поменяем a на b, p на q и x на y, так что теперь мы всегда ищем точку x ниже точки перегиба<xs>и на выпуклой кривой. Первоначальная оценка для x производится с:

который оказывается ниже истинного значения для x:Таким образом, итерация Ньютонабудет плавно сходиться на x без проблем, вызванных превышением и т.д.

Если a и b больше 1, но a+b слишком малы, чтобы использовать другие методы, упомянутые выше, мы приступаем к следующему. Обратите внимание, что в неполной бета-версии есть точка перегиба<xs =(1-a)/(2-a -b)>. Поэтому, если<p > Ix(a,b)>мы поменяем a на b, p на q и x на y, так что теперь мы всегда ищем точку x ниже точки перегиба<xs>и на вогнутой кривой. Первоначальная оценка для x производится с:

которые можно несколько улучшить до:

когда b и x малы (я использовал b

Implementation of inverses on the a and b parameters

Эти четыре функции имеют общую реализацию.

Сначала вычисляется начальное приближение дляaилиb: там, где это возможно, используется расширение Корниш-Фишера для отрицательного биномиального распределения, чтобы попасть в пределах около 1 результата. Однако, когдаaилиbочень малы, расширение Корниша Фишера не может быть использовано, в этом случае первоначальное приближение выбрано так, что I_x(a, b) находится вблизи середины диапазона [0,1].

Эта исходная догадка затем используется в качестве начального значения для общего алгоритма поиска корней. Алгоритм быстро сходится на корне после того, как он был скоблен, но для скобки корня может потребоваться несколько итераций. Лучшее первоначальное приближение дляaилиbзначительно улучшило бы эти функции: в настоящее время для поиска корня требуется 10-20 неполных вызовов бета-функции.