Практические вычисления интегрального F - распределения

[Andrey_Kireew]

Использование статистических таблиц в компьютерной программе не всегда удобно, особенно, в случае F - распределения, зависящего сразу от 2-х параметров:

$$F(x | k,m)=\frac{\Gamma ((k+m)/2)}{\Gamma (k/2)\Gamma (m/2)}\int\limits_{0}^{x}\frac{(k/m)^{k/2}f{(k-1)/2}}{(1+t\cdot k/m)^{(k+m)/2}}dt$$

Естественно возникает вопрос: есть ли для этой функции какие либо аппроксимации или разложения, удобные для практических вычислений?

Для других распределений (нормального, Хи - квадрат) такие формулы представлены в большом разнообразии, и найти их не составляет особого труда. Но по данному вопросу информации практически нет. Единственное, что удалось найти - книга Абрамовица М., Стиган И. Справочник по специальным математическим функциям. - М. "Наука", 1979г. В ней представлено несколько сомнительных формул. Но к ним нет никаких описаний. Главное - непонятно какова их точность.

Буду благодарен любой полезной информации по данному вопросу

[novichok2018]

Из под интеграла действительно можно вынести всё что сверху?
И в адрес такой уважаемой профессионалами книги, как Абрамовиц-Стиган, я бы не стал высказываться, что там сомнительные формулы. Что теперь есть более новые и полные справочники - не влияет.

[arseniiv]

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

-- Сб дек 07, 2019 19:29:48 --

Можно заметить, что всего этого я раньше не знал и просто нашёл чуть-чуть поискав.

-- Сб дек 07, 2019 19:31:28 --

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

[Andrey_Kireew]

arseniiv в сообщении #1429185 писал(а):

Можно заметить, что всего этого я раньше не знал и просто нашёл чуть-чуть поискав

оно и заметно ...

ничего там не выражается, бета-распределение вычислить так же сложно, как и F - распределение, интеграл от такой замены никуда не исчезает. Напрямую, через гамма-функцию и численное интегрирование, вычислять слишком громоздко. А никаких разложений или аппроксимаций для бета-распределения, разумеется, я точно также не нашел. Потому и написал, подумав, что может кто знает решение проще. И у меня есть подозрение, что решение это существует

[arseniiv]

Andrey_Kireew в сообщении #1429191 писал(а):

Напрямую, через гамма-функцию и численное интегрирование, вычислять слишком громоздко.

Так на практике гамма-функция (обычная ли, логарифм ли её) и вычисляется с некоторыми попущениями, чтобы побыстрее. И интегрирование в явном виде не понадобится, только гамма.

-- Сб дек 07, 2019 20:18:37 --

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

[Andrey_Kireew]

Да потому, что пересмотрел уже всё 100 раз ... нет там ничего хорошего. Самое близкое к тому что нужно, как я уже писал, есть в книжке Абрамовица М., Стиган И ... (к стати, бета-распределения там не обнаружил), но там не указана точность этих приближенных формул, как это обычно делается. Без этого их использовать проблематично.

На счёт вычисления гамма-функции вопросов нет, с этим всё предельно ясно. А как вычислять интегральную бета-функцию, точно так же, нигде не сказано. И как же тогда не понадобится интегрирование в явном виде?

Я к тому, что сама по себе, эта замена мало что даёт, так как нужен способ вычисления бета-распределения.
По любезно предоставленной Вами ссылке я нашел изящное решение

Цитата:

An incomplete beta value can be calculated as: Value = BetaDist(x, a, b) * Exp(GammaLn(a) + GammaLn(b) − GammaLn(a + b)).
These result follow from the properties listed above.

но в math.h нет никакой BetaDist(), хотя гамма-функция там и есть.

Получается, вопрос как был, так и остался. Кроме численного интегрирования я пока вариантов не вижу

[arseniiv]

BetaDist определяется строчкой выше.

Кстати если пишете на C++ и компилятор поддерживает C++17, есть и бета-функция (не неполная). Если поискать по c math library beta, можно найти код для C аж с объяснением того, как он получен: https://codeplea.com/incomplete-beta-function-c (найти внезапно, потому что предполагалось найти код из какой-нибудь целой библиотеки, а не отдельный пост; и например код из Boost находится, но он не для C, а для плюсов).

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

[Andrey_Kireew]

arseniiv в сообщении #1429196 писал(а):

можно найти код для C аж с объяснением того, как он получен: https://codeplea.com/incomplete-beta-function-c

Спасибо arseniiv, это похоже на то, что нужно. Буду переводить, я по английски читаю плохо.

[novichok2018]

Про формулу ещё раз - там действительно интеграл по t, но в числителе под интегралом и в самой функции ничего от $t$ не зависит?
Про F распределение. Есть очень хороший современный справочник по теме в 2 томах.
Джонсон и др. Одномерные непрерывные распределения. 2012.
Во 2 томе в пар. 27 всё про F распределения.

[Andrey_Kireew]

novichok2018 в сообщении #1429283 писал(а):

Про формулу ещё раз - там действительно интеграл по t, но в числителе под интегралом и в самой функции ничего от $t$ не зависит?

Просто в формуле ошибка, вместо f должно быть t (иначе возникает вопрос, что же такое f), я потом только заметил, исправлять было уже поздно.
Спасибо за ссылку, правда скачать пока её не получилось.

Хотелось бы отдельно поблагодарить arseniiv, за предоставленную ссылку. Алгоритм Lentz's оказался на удивление хорош. Я его тут "погонял" немного на разных числовых значениях, получается максимальная относительная погрешность $10^{-7}$, а средняя $10^{-12}$. Если уменьшить критерий останова с $10^{-8}$ до $10^{-9}$, точность возрастает в 10 раз. Потом она уже не увеличивается, видимо, сказывается точность вычисления гамма-функции. Работает быстро, обычно сходится за 4-5 итераций. В редких случаях (примерно один из 20), количество итераций возрастает до нескольких десятков. Максимум - это 50.

К стати, гамма функция реализована только в С++11, у меня её к сожалению тоже не было. Так, что исходный код пришлось немного дополнить.
Но это всё мелочи