Оценить определенный интервал нормального распределения

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

Однако мне интересно, есть ли хорошая формула для . Или каково «современное» приближение для этой проблемы.$N(c_{-} \leq x < c_{+}| \mu, \sigma^2)$

Это зависит от того, что именно вы ищете . Ниже приведены краткие сведения и ссылки.

Большая часть литературы по приближениям сосредоточена вокруг функции

$$Q(x) = \int_x^\infty \frac{1}{\sqrt{2\pi}} e^{-\frac{u^2}{2}} \, \mathrm{d}u$$

для . Это связано с тем, что предоставленная вами функция может быть разложена как простая разница вышеупомянутой функции (возможно, скорректированная на константу). Эта функция называется многими именами, включая «верхний хвост нормального распределения», «правый нормальный интеграл» и « функцию Гаусса », чтобы назвать несколько. Вы также увидите приближения к отношению Миллса , которое равно где - это гауссовский pdf.$x > 0$$Q$ φ(х)=(2π)-1/2е-х2/

$$R(x) = \frac{Q(x)}{\varphi(x)}$$ $\varphi(x) = (2\pi)^{-1/2} e^{-x^2 / 2}$

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

вычислительный

Де-факто стандартом для вычисления функции или связанной с ней дополнительной функции ошибки является$Q$

WJ Cody, Рациональные чебышевские приближения для функции ошибок , Math. Комп. , 1969, с. 631-637.

Каждая (уважающая себя) реализация использует этот документ. (MATLAB, R и т. Д.)

«Простые» приближения

Абрамовиц и Стегун имеют один, основанный на полиномиальном разложении преобразования входных данных. Некоторые люди используют это как "высокоточное" приближение. Мне не нравится это для этой цели, так как это ведет себя плохо около нуля. Например, их приближение не дает , что, я думаю, является большим нет-нет. Иногда из-за этого случаются плохие вещи .$\hat{Q}(0) = 1/2$

Боржессон и Сундберг дают простое приближение, которое работает довольно хорошо для большинства приложений, где требуется лишь несколько цифр точности. Абсолютная относительная погрешность никогда не хуже , чем на 1%, что вполне неплохо , учитывая его простоту. Основное приближение: и их предпочтительный выбор констант: и . Эта ссылкаa=0,339b=

$$\hat{Q}(x) = \frac{1}{(1-a) x + a \sqrt{x^2 + b}} \varphi(x)$$ $a = 0.339$$b = 5.51$

П.О. Боржессон и С.Е. Сундберг. Простые аппроксимации функции ошибок Q (x) для приложений связи . IEEE Trans. Commun. COM-27 (3): 639–643, март 1979 г.

Вот график его абсолютной относительной погрешности.

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

Вы также можете посмотреть на

В. Брайк. Равномерное приближение к правому нормальному интегралу . Прикладная математика и вычисления , 127 (2-3): 365–374, апрель 2002.

Продолжение Лапласа

У Лапласа есть красивая непрерывная дробь, которая выдает последовательные верхнюю и нижнюю границы для каждого значения . Это, с точки зрения отношения Миллса,$x > 0$

$$R(x) = \frac{1}{x+}\frac{1}{x+}\frac{2}{x+}\frac{3}{x+}\cdots ,$$

где обозначение, которое я использовал, является довольно стандартным для непрерывной дроби , то есть . Это выражение не очень быстро сходится при малых значениях и расходится при .x x $1/(x+1/(x+2/(x+3/(x+\cdots))))$$x$$x = 0$

Эта непрерывная дробь фактически дает множество «простых» границ которые были «переоткрыты» в середине-конце 1900-х годов. Легко видеть, что для непрерывной дроби в «стандартной» форме (т. Е. Составленной из положительных целых коэффициентов), усечение дроби в нечетных (четных) терминах дает верхнюю (нижнюю) границу.$Q(x)$

Следовательно, Лаплас немедленно сообщает нам, что обе из которых являются «переоткрытыми» в середине 1900 -х годов. С точки зрения функции это эквивалентно Альтернативное доказательство этого с помощью простого интегрирования по частям можно найти в работе С. Резника, « Приключения в случайных процессах» , Birkhauser, 1992, в главе 6 (Броуновское движение). Абсолютная относительная погрешность этих границ не хуже, чем , как показано в этом связанном ответе .Q

$$\frac{x}{x^2 + 1} < R(x) < \frac{1}{x} \>,$$ $Q$ $$\frac{x}{x^2 + 1} \varphi(x) < Q(x) < \frac{1}{x} \varphi(x) .$$ $x^{-2}$

В частности, обратите внимание, что из приведенных выше неравенств немедленно следует, что . Этот факт может быть установлен с использованием правила Лопиталя. Это также помогает объяснить выбор функциональной формы приближения Боржессона-Сундберга. Любой выбор поддерживает асимптотическую эквивалентность как . Параметр служит «коррекцией непрерывности» вблизи нуля.$Q(x) \sim \varphi(x)/x$$a \in [0,1]$$x \to \infty$$b$

Вот график функции и двух оценок Лапласа.$Q$

CI. У К. Ли есть статья начала 1990-х годов, в которой делается «коррекция» для малых значений . Видеть$x$

CI. К. Ли. На Лапласе продолжается дробь для нормального интеграла . Энн. Текущий месяц Statist. Математика , 44 (1): 107–120, март 1992.

---

Вероятность Дарретта : теория и примеры дают классические верхние и нижние оценки на страницах 6–7 3-го издания. Они предназначены для больших значений (скажем, ) и асимптотически тесны.$Q(x)$$x$$x > 3$

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

Полагаю, я опоздал, герой, но я хотел прокомментировать сообщение кардинала, и этот комментарий стал слишком большим для предполагаемой рамки.

Для этого ответа я предполагаю, что ; соответствующие формулы отражения могут быть использованы для отрицательного х .$x >0$$x$

Я более привык иметь дело с функцией ошибок сам, но я попытаюсь пересмотреть то, что я знаю, с точки зрения отношения Миллса R ( x ) (как определено в ответе кардинала).$\mathrm{erf}(x)$$R(x)$

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

$|x|$

$$R(x)=\sqrt{\frac{\pi}{2}}\exp\left(\frac{x^2}{2}\right)-x\sum_{j=0}^\infty\frac{2^j j!}{(2j+1)!}x^{2j}$$

$x^{2j}$$c_j=\frac{2^j j!}{(2j+1)!}$$c_0=1$$c_{j+1}=\frac{c_j}{2j+3}$

---

$|x|$

Ленц , Томпсон и Барнетт вывели алгоритм численной оценки непрерывной дроби как бесконечного произведения, который более эффективен, чем обычный подход вычисления непрерывной дроби «назад». Вместо отображения общего алгоритма я покажу, как он специализируется на вычислении отношения Миллса:

$\displaystyle Y_0=x,\,C_0=Y_0,\,D_0=0$
$\text{repeat for }j=1,2,\dots$

$$D_j=\frac1{x+jD_{j-1}}$$ $$C_j=x+\frac{j}{C_{j-1}}$$ $$H_j=C_j D_j$$ $$Y_j=H_j Y_{j-1}$$ $\text{until }|H_j-1| < \text{tol}$
$\displaystyle R(x)=\frac1{Y_j}$

$\text{tol}$

CF полезен, когда упомянутые ранее серии начинают медленно сходиться; вам придется поэкспериментировать с определением подходящей «точки останова» для перехода от серии к CF в вашей вычислительной среде. Существует также альтернатива использования асимптотического ряда вместо CF Лапласа, но мой опыт показывает, что CF Лапласа достаточно хорош для большинства приложений.

---

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

$$R(x)\approx \frac{\sqrt{2\pi}+x(\pi-2)}{2+x\sqrt{2\pi}+x^2(\pi-2)}$$

$1.84\times 10^{-2}$$x$

(Этот ответ первоначально появился в ответ на аналогичный вопрос, впоследствии закрытый как дубликат. ФП хотел только «реализацию» гауссовского интеграла, а не обязательно «современное состояние». В его комментариях стало очевидно, что относительно простой короткая реализация будет предпочтительнее.)

---

$-8.5$$+8.5$

Версия MatLab (с соответствующими атрибутами) доступна по адресу http://people.sc.fsu.edu/~jburkardt/m_src/asa005/alnorm.m . Полностью недокументированная версия исходного кода на Фортране появляется на сайте «Поиск кода Koders» (sic).

Много лет назад я перенес это в AWK. Эта версия может быть более подходящей для современного разработчика для портирования из-за ее C-подобного (а не Fortran) синтаксиса и некоторых дополнительных комментариев, которые я вставил при разработке и тестировании, потому что мне нужно было повысить ее точность. Похоже, ниже.

Для тех, у кого нет большого опыта переноса научного / математического / статистического кода, несколько советов : одна типографская ошибка может привести к серьезным ошибкам, которые трудно обнаружить. (Поверьте мне, я сделал их много.) Всегда, всегда создавайте тщательный и исчерпывающий тест. Поскольку нормальная функция интеграла / гауссова интеграция / ошибка доступна во многих таблицах и большом количестве программного обеспечения, просто и быстро составить таблицу огромного числа значений вашей перенесенной функции и систематически сравнивать (т. Е. С компьютером, а не на глаз) значения для исправления. Вы можете увидеть такой тест в начале моего кода: он создает таблицу значений в -8,5: 8,5 (на 0,1), которую можно передать (через STDOUT) в другую программу для систематической проверки.

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

$0$$-\infty$$x$$\mu$$\sigma$$z = (x-\mu)/\sigma$alnorm

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

alnorm$1 - \Phi(z)$$z\ge 1$alnorm

$4 \times 10^{-11}$ $z=16$$z$$z=\sqrt(2\times 708) \approx 37.6$

alnorm[-6.0]$9.865\ 876\ 450\ 315E-10$$\frac{1}{2}\text{erfc}(3\sqrt{2})$$9.865\ 876\ 450\ 377E-10$

UPPER_TAIL_IS_ZERO15.16.$Z$$15$$16$

#----------------------------------------------------------------------#
#   ALNORM.AWK
#   Compute values of the cumulative normal probability function.
#   From G. Dallal's STAT-SAK (Fortran code).
#   Additional precision using asymptotic expression added 7/8/92.
#----------------------------------------------------------------------#
BEGIN {
    for (i=-85; i<=85; i++) {
        x = i/10
        p = alnorm(x, 0)
        printf("%3.1f %12.10f\n", x, p)
    }
    exit
}
function alnorm(z,up,    y,aln,w) {
#
#    ALGORITHM AS 66 APPL. STATIST. (1973) VOL.22, NO.3:
#    Hill,  I.D.  (1973).  Algorithm AS 66.  The normal  integral.
#                          Appl. Statist.,22,424-427.
#
#    Evaluates the tail area of the standard normal curve from
#    z to infinity if up, or from -infinity to z if not up.
#
#    LOWER_TAIL_IS_ONE, UPPER_TAIL_IS_ZERO, and EXP_MIN_ARG
#    must be set to suit this computer and compiler.

    LOWER_TAIL_IS_ONE = 8.5     # I.e., alnorm(8.5,0) = .999999999999+
    UPPER_TAIL_IS_ZERO = 16.0   # Changes to power series expression
    FORMULA_BREAK = 1.28        # Changes cont. fraction coefficients
    EXP_MIN_ARG = -708          # I.e., exp(-708) is essentially true 0

    if (z < 0.0) {
        up = !up
        z = -z
    }
    if ((z <= LOWER_TAIL_IS_ONE) || (up && z <= UPPER_TAIL_IS_ZERO)) {
        y = 0.5 * z * z
        if (z > FORMULA_BREAK) {
            if (-y > EXP_MIN_ARG) {
                aln = .398942280385 * exp(-y) / \
                  (z - 3.8052E-8 + 1.00000615302 / \
                  (z + 3.98064794E-4 + 1.98615381364 / \
                  (z - 0.151679116635 + 5.29330324926 / \
                  (z + 4.8385912808 - 15.1508972451 / \
                  (z + 0.742380924027 + 30.789933034 / \
                  (z + 3.99019417011))))))
            } else {
                aln = 0.0
            }
        } else {
            aln = 0.5 - z * (0.398942280444 - 0.399903438504 * y / \
              (y + 5.75885480458 - 29.8213557808 / \
              (y + 2.62433121679 + 48.6959930692 / \
              (y + 5.92885724438))))
        }
    } else {
        if (up) {   # 7/8/92
            # Uses asymptotic expansion for exp(-z*z/2)/alnorm(z)
            # Agrees with continued fraction to 11 s.f. when z >= 15
            # and coefficients through 706 are used.
            y = -0.5*z*z
            if (y > EXP_MIN_ARG) {
                w = -0.5/y  # 1/z^2
                aln = 0.3989422804014327*exp(y)/ \
                    (z*(1 + w*(1 + w*(-2 + w*(10 + w*(-74 + w*706))))))
                # Next coefficients would be -8162, 110410
            } else {
                aln = 0.0
            }
        } else {
            aln = 0.0
        }
    }
    return up ? aln : 1.0 - aln
}
### end of file ###