Поиск хи-квадрат по квантилю

[mitrik]

Доброго времени суток! Имеется следующая задача : на каждой итерации есть число, вещественное, оно является 5% квантилем для одного хи-квадрат и 95% квантилем для другого хи-квадрат. Эти распределения имеют одинаковое количество степеней свободы, это число известно, отличаются только смещением, надо найти разницу между их медианами. Также очень интересует вопрос сходимости этой самой разности

[Евгений Машеров]

Вообще-то единственный параметр распределения хи-квадрат это число степеней свободы. Так что если это два разных распределения с одинаковым параметром "число степеней свободы", то хотя бы одно из них не хи-квадрат.

[mitrik]

Хи-квадраты нецентральные, поэтому отличаются параметрами смещения

[GAA]

mitrik, пожалуйста, приведите значение параметра число степеней свободы.

[mitrik]

GAA, этот параметр известен, но это также переменная, равная текущей итерации, поэтому числового значения я, к сожалению, представить не могу

[GAA]

mitrik, в чем конкретно состоит задача?

Если нужно для некоторой заданной последовательности чисел найти «разницу» между медианами, то в этом случае в Сети можно отыскать калькулятор и воспользоваться им, например Keisan online calculator.

Если нужно написать программу, то
1. Нецентральное $\chi^2$ распределение есть в MatLab. Можно попробовать написать скрипт в этой системе. [В справке есть ссылка на английское издание Абрамовиц М., Стиган И. Справочник по специальным функциям. Может и в русском издании есть нужные сведения, см. раздел 26.4 "Распределение хи-квадрат" в издании 1979 г.]
2. В разделе «Назначение таблиц распределения $\chi^2$ и примеры их использования» книги Большев Л.Н., Смиронов Н.В. Таблицы математической статистики, 1983 приведены аппроксимации, в частности Пирсона, и ссылки. Можно написать программу для быстрого грубого вычисления.

[Евгений Машеров]

Я не вижу иного пути, как вычислить значения двух параметров нецентральности, затем найти для этих распределений медианы и получить их разность.
Работ по оценке параметров нецентрального $\chi^2$ немало, но в основном они посвящены оценке по выборке. Оценка по значениям квантилей мне не попадалась. В принципе, функция распределения известна и можно считать по ней - но она страшненькая, бесконечная сумма, а под знаком суммы интеграл живёт. Аппроксимаций, удобных для оценивания, для обратной функции я не нашёл, возможно, и есть такие, но не увидел.
В порядке грубой оценки - берём аппроксимацию для функции распределения нецентрального $\chi^2$, в основном они относятся к двум типам - через нормальное распределение и через центральное $\chi^2$, первое мне представляется удобнее для данной цели, зная, что это 5% и 95% квантили, получаем значение аргумента аппроксимирующей нормальной функции, -2 и +2 соответственно (я думаю, что - 1.96 и + 1.96 при такой грубой аппроксимации это бессмысленное пижонство), и затем через аргумент и его выражение через параметры распределения нецентрально $\chi^2$ получаем эти параметры.
Вот есть (очень грубая, полученная приравниванием моментов) аппроксимация
$F(x;\nu,\lambda) =\Phi(\frac{x -\nu -\lambda} {\sqrt{2(\nu+2\lambda)}})$
Приравнивая аргумент к плюс или минус двум и решая, зная число степеней свободы и значение квантиля, относительно параметра нецентральности, находим его. Затем той же аппроксимацией или более точной оцениваем медиану и находим разность (для данной аппроксимации разность медиан будет попросту разностью параметров нецентральности, в более точных используются нелинейные преобразования. Эту аппроксимацию нашёл в Continuous univariate distributions. Vol.2 Johnson N.L., Kotz S., Balakrishnan N. N.-Y.: Wiley, 1995, там есть и другие.

[Евгений Машеров]

У Абрамовица и Стигана приведена аппроксимация обратной функции нецентрального $\chi^2$-распределения
Если $Q(\chi^2_\lambda|\nu,\lambda)=p$ и $Q(\chi^2|\nu^*)=p$
где $a=\nu+\lambda$, $b=\frac \lambda {\nu+\lambda}$, $\nu^*=\frac a {1+b}=\frac {(\nu+\lambda)^2}{\nu+2\lambda}$
то $\chi^2_\lambda=(1+b)\chi^2$
К сожалению, не указана погрешность и границы применимости аппроксимации. Впрочем, можно надеяться, что ошибки аппроксимации будут одного знака и порядка и отчасти компенсируют друг друга после вычитания.
Так что можно, получив значения 5% и 95% квантилей для нецентральных распределений, найти для них значения параметров нецентральности (увы, число степеней свободы в аппроксимирующем хи-квадрат зависит от нецентральности аппроксимируемого, так что решается нелинейное уравнение с неизвестным в качестве аргумента неэлементарной функции), затем подставить в выражение для 50% квантилей нецентрального и найти искомую разность.

[mitrik]

Для поиска параметров нецентральности решил воспользоваться аппроксимацией Фишера $\chi_p(f) = 0.5(u_p + \sqrt{2f-1})^2$, $u_p $ - искомый квантиль для нормального распределения. Для центрального хи-квадрат и перехода к нему от нецентрального использовалась формула Патнайка ($\chi(f') = \frac{f + \lambda}{f + 2\lambda}\chi(f ,\lambda)$), где ($f' = \frac{(f + 9)^2}{f + 2\lambda}$). Однако с ростом числа степеней свободы и числа данного в условии, на месте параметра нецентральности появляются комплексные числа, с чем это может быть связано?

Уравнение решалось относительно параметра нецентральности в WolphramAlpha

Уравнение имеет вид : $0.5(u_p + \sqrt{2\frac{(f + 9)^2}{f + 2\lambda} - 1})^2 = p \frac{f + \lambda}{f + 2\lambda}$, где $p$ - данное в условии число

[Евгений Машеров]

Похоже, Вы вышли за границы применимости этой аппроксимации. "Аппроксимации аппроксимации", я бы сказал, потому как Патнайк приближает нецентральный $\chi^2$ центральным, а Фишер центральный $\chi^2$ нормальным. И где-то было столь тонко, что порвалось.
С другой стороны - смущает меня выражение для приближения Патнайка. Эннеада там откуда взялась, в смысле девятка? Согласно Кендаллу и Стьюарту ("Статистические выводы и связи", М.: Наука, 1973, с. 307, формулы 24.21) степени свободы для приближающего нецентральное центрального будут (в Ваших обозначениях) не $f' = \frac{(f + 9)^2}{f + 2\lambda}$, а вовсе даже $f' = \frac{(f + \lambda)^2}{f + 2\lambda}$

[mitrik]

Использовал книгу Кобзаря А.И. ("Прикладная математическая статистика. Для инженеров и научных работников", М.:ФИЗМАТЛИТ, 2006, с. 81, формула 99). Не использовал описанную Вами аппрокисмацию по Кендаллу и Стюарту из-за излишней громоздкости(хотя тут это вообще нормальное явление) и множественности корней, так как задача требует быстроты и хорошей работы в динамике

[Евгений Машеров]

У меня по поводу названной книги сложилось не вполне благоприятное впечатление. Там слишком много опечаток, а так как это справочник, а не оригинальное исследование и не учебник, где можно отследить логику рассуждения и восстановить правильный вид, то выявить ошибочность нельзя. Разве что интуитивно - что это ещё за магическое число 9?
В своё время подобного рода справочники тщательно вычитывались мощным коллективом редакторов и корректоров, но тут, кажется, редактор один, он же и составитель, а на корректорах издательства изрядно экономят.
Боюсь, что пользоваться ею неосторожно, и максимум, на что она годится - подыскав с её помощью нужную формулу, найти, откуда она взялась, отыскать ссылку на первоисточник и там найти правильный вид.
Что тамошнее (Кобзаря) приближение неверное, впрочем, можно заметить не только из эстетических соображений, нецентральное распределение совпадает с центральным, если параметр нецентральности $\lambda=0$. В этом случае "вырожденной аппроксимации" $f'=f$, а никак не, как даёт Кобзарь, $f'=f+18+\frac {81}f$

[mitrik]

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

1) В вычислении с помощью

Евгений Машеров в сообщении #1089272 писал(а):

$F(x;\nu,\lambda) =\Phi(\frac{x -\nu -\lambda} {\sqrt{2(\nu+2\lambda)}})$

не очень понимаю принцип и порядок действий, так как функция Лапласа ограничена, соответственно и количество решений будет большим

2) Аппроксимации обратной функции не нравятся тем, что придется решать нелинейное уравнение с невеселым интегралом, что потребует численных методов, а в задаче важно быстродействие

3) Метод аппроксимации через Патнайка ведет к уравнению 4ой степени, которое даже WolphramAlpha иногда отказывается решать

Посему я в глубоком раздумьи : искать ли мне улучшение текущих методов или плюнуть и реализовать какой-либо из них

P.S. реализация будет на С

[Евгений Машеров]

Аппроксимации, основанные на нормальном распределении, все асимптотические, и при малом $\nu$ могут быть очень плохими (в частности, потому, что при больших $\nu$ мы вряд ли попадём туда, где левый хвост отрицателен, чего у $\chi^2$ быть не может. А так как у Вас на первых шагах число степеней свободы мало - то свою рекомендацию снимаю.
Решать через интеграл - да, медленно.
Я бы пользовался Патнайком. Только там не четвёртая степень, там хуже. Там фигурирует центральная $\chi^2$ с числом степеней свободы, зависящим от параметра нецентральности аппроксимируемого. То есть даже не алгебраическое уравнение четвёртой степени, а нелинейное общего вида. Решать его надо численно, делением пополам (бисекциями), regula falsi (интерполяцией) или Ньютоном. Причём надо уметь вычислять квантили $\chi^2$ для разных степеней свободы.
Вот некоторые алгоритмы, но найти я их нашёл, но не изучил:
фортрановский (достаточно несложно переписываемый на С, и даже некогда, помнится, я его использовал, но исходники утерял)

   FUNCTION CHISQD(P, N)

c*********************************************************************72

      DIMENSION C(21), A(19)
      DATA C(1)/1.565326E-3/, C(2)/1.060438E-3/,
     &  C(3)/-6.950356E-3/, C(4)/-1.323293E-2/,
     &  C(5)/2.277679E-2/, C(6)/-8.986007E-3/,
     &  C(7)/-1.513904E-2/, C(8)/2.530010E-3/,
     &  C(9)/-1.450117E-3/, C(10)/5.169654E-3/,
     &  C(11)/-1.153761E-2/, C(12)/1.128186E-2/,
     &  C(13)/2.607083E-2/, C(14)/-0.2237368/,
     &  C(15)/9.780499E-5/, C(16)/-8.426812E-4/,
     &  C(17)/3.125580E-3/, C(18)/-8.553069E-3/,
     &  C(19)/1.348028E-4/, C(20)/0.4713941/, C(21)/1.0000886/
      DATA A(1)/1.264616E-2/, A(2)/-1.425296E-2/,
     &  A(3)/1.400483E-2/, A(4)/-5.886090E-3/,
     &  A(5)/-1.091214E-2/, A(6)/-2.304527E-2/,
     &  A(7)/3.135411E-3/, A(8)/-2.728484E-4/,
     &  A(9)/-9.699681E-3/, A(10)/1.316872E-2/,
     &  A(11)/2.618914E-2/, A(12)/-0.2222222/,
     &  A(13)/5.406674E-5/, A(14)/3.483789E-5/,
     &  A(15)/-7.274761E-4/, A(16)/3.292181E-3/,
     &  A(17)/-8.729713E-3/, A(18)/0.4714045/, A(19)/1./

      IF (N-2) 10, 20, 30
10    CHISQD = GAUSSD(.5*P)
      CHISQD = CHISQD*CHISQD
      RETURN
20    CHISQD = -2. * ALOG(P)
      RETURN
30    F = N
      F1 = 1. / F
      T = GAUSSD(1.-P)
      F2 = SQRT(F1) * T
      IF ( N.GE.(2+INT(4.*ABS(T)))) GO TO 40
      CHISQD=(((((((C(1)*F2+C(2))*F2+C(3))*F2+C(4))*F2
     &  +C(5))*F2+C(6))*F2+C(7))*F1+((((((C(8)+C(9)*F2)*F2
     &  +C(10))*F2+C(11))*F2+C(12))*F2+C(13))*F2+C(14)))*F1 +
     &  (((((C(15)*F2+C(16))*F2+C(17))*F2+C(18))*F2
     &  +C(19))*F2+C(20))*F2+C(21)
      GO TO 50
40    CHISQD = (((A(1)+A(2)*F2)*F1+(((A(3)+A(4)*F2)*F2
     &  +A(5))*F2+A(6)))*F1+(((((A(7)+A(8)*F2)*F2+A(9))*F2
     &  +A(10))*F2+A(11))*F2+A(12)))*F1 + (((((A(13)*F2
     &  +A(14))*F2+A(15))*F2+A(16))*F2+A(17))*F2*F2
     &  +A(18))*F2+A(19)
50    CHISQD = CHISQD*CHISQD*CHISQD*F

      RETURN
      END

http://people.sc.fsu.edu/~jburkardt/f77 ... ms451.html

и ещё один, более сложный (но, как уверяют авторы, более точный) и тоже на фортране
https://wiki.hepg.sdu.edu.cn/doc/cernli ... antile.pdf

Кроме того, есть вариант, имея на руках таблицы распределения, выписать оттуда значения для разных степеней свободы для 5%, 50% (медианы) и 95%, и брать оттуда (интерполируя, если в аппроксимации появится дробное $\nu$). Возможно, я так бы и поступил.

[mitrik]

Число степеней свободы, при котором начнет выполнятся данный блок, как правило превышает 10-15, при таком количестве уже можно пользоваться аппроксимацией через нормальное, где, если я правильно понимаю
$F(x;\nu,\lambda) =\Phi(\frac{x -\nu -\lambda} {\sqrt{2(\nu+2\lambda)}})$ эту формулу $x_k -$ доступное мне на $k-$ой итерации число, а приравнивать аргумент функции Лапласа нужно к плюс или минус двойке?