Вычисление функции интеграла вероятности и ей обратной Текст научной статьи по специальности «Математика»

УДК 518.5

Г.С. ТЕСЛЕР, ЗЫ ХАК ЗУНГ

ВЫЧИСЛЕНИЕ ФУНКЦИИ ИНТЕГРАЛА ВЕРОЯТНОСТИ И ЕЙ ОБРАТНОЙ

Abstract: In the article therehase been proposed the algorithms for effective calculating probabilities integral function and the reverse bit founded on decomposing these functions in row of closure errors on polynomials Erpith. For calculating probabilities integral quantiles there has been received the functions decomposing in row of closure errors on the basis of converting the probabilities integral functions uhat requires considerably lesser expenses than Jailor’s decomposing in raw/ On the basis of this decomposing there have been received the iterative formulae of cowergency arbitrary arrangement. Special attention has been spared to initial estimates on the basis of segment approximation? To the problems of tabulating the functions and computing mutually-inverse functions. Shere have been represented the materials of modelling the proposed algorithms on a computer.

Key words: probabilities integral, quantile, approximation, modelling, error.

Анотація: В статті запропоновані алгоритми для ефективного обчислення функції інтеграла імовірностей та її оберненої, які базуються на розкладі цих функцій у ряд за нев’язками по поліномах Ерміта. Для обчислення квантилей інтеграла імовірностей отримано розклад функцій у ряд нев’язок на основі обернення розкладу функції інтеграла імовірностей, що потребує значно менших витрат, ніж розклад у ряд Тейлора. На базі цього розкладу отримані ітераційні формули довільного порядку збіжності. Особливу увагу присвячено початковим наближенням на основі сегментної апроксимації, питанням табулювання функцій та обчислення взаємно-обернених функцій. Представлені матеріали моделювання запропонованих алгоритмів на комп’ютері.

Ключові слова: інтеграл імовірностей, квантиль, апроксимація, моделювання, похибки.

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

Ключевые слова: интеграл вероятностей, квантиль, аппроксимация, моделирование, погрешность.

1. Введение

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

Это связано с тем, что данные функции непосредственно связаны с нормальным (гауссовским) распределением случайной величины, которые во многих случаях являются адекватной моделью многих процессов, происходящих в живой и искусственной природе. Они широко используются в теории вероятности, математической статистике, теоретической и математической физике, описывающей процессы теплопроводности, диффузии, броуновского движения, кинетической теории газов и т.д., в теории связи (при исследовании влияния шума на полезный сигнал), при решении дифференциальных и интегральных уравнений методами Монте Карло, в теории чисел (распределение значений “коротких” чисел Вейля и других разделах теоретической и прикладной математики, а также в теории надёжности и многих других приложениях).

Описание функции интеграла вероятностей либо их таблиц включено в большинство известных справочников и монографий, посвящённых описанию специальных математических функций [1 - б]. Имеется обширная литература с подробными таблицами данных функций [7, В]. В [В]

приведена обширная библиография по этой тематике. Имеются также отдельные монографии, посвящённые функции интеграла вероятности [9, 26], и множество таблиц.

В работах [10 - 13] приведены различные методы и алгоритмы аппроксимации функции вероятности, а также обширная библиография.

Особо хочется отметить тот факт, что в шестидесятые годы прошлого столетия в работах [14 -15] ввели в рассмотрение функцию у = inverf (х), являющуюся обратной по отношению функции интеграла вероятности. Это не означает, что ранее в неявном виде данные функции использовались под названием квантиля распределения и процентных точек распределения.

Для эффективного вычисления как прямых, так и обратных функций на ЭВМ, Г.С. Теслером было предложено использовать адаптивные аппроксимации, основанные на разложении функций по невязкам (см. справочник [13] и приведенную в нём библиографию).

В настоящей статье предлагается дальнейшее развитие этих исследований для функций

2 х 1 X

erf (х) = —1= [e- dt , Ф0(х) = . [e-t ^dt и им обратным. Отметим, что erf (х ) = 2Ф 0(хТ2) .

4ж о л/2р 0

В свою очередь функция

1 х 1 1

Ф(х) = -= je-t2/2dt = -[1 - erf (х/ V2)] = - + Фо (х).

V2p 2 2

2. Вычисление erf (х) и Ф 0( х)

Разложение erf (х + у) в ряд Тейлора имеет вид

¥ 1 к к ¥ 1 к+1 к+1

erf (х+у) = X -vjerf ( у^ у=х • -г: = erf (х) + £ тз+тerf ( у^ у=х •

к=0

Учтём, что производная интеграла вероятностей

dyk " ^'|у=х к! - ' к=о dyk+1 ^ ~ '|у=х (к +1)!

. derf (х) 2 - х2

j( х) = ^ J=~j= e х ;

dх

dn+1 2 2

jn)( х) = — erf (х) = (-1) n-= Нп (х)*-х = (-1) nHn (х). j х), n = 0,1,2...;

dхn+1 yjp

(pn+1 (х) + хфп (х) + n (pn-1 (х) = 0, j(0) (х) = j( х),

Hn (х) - многочлены Эрмита.

Отметим, что

Hn (х) = (-1) ne* —e -х dх

Произведя несложные преобразования и учтя, что невязка z0 = inverf (у0) - х = х0 - х и у0 приближённое значение y = erf (х) на заданном интервале, получим

erf ( х) = erf (inverf (у0) - ^0) = у0 —т= e4mve,f(у0)) .£

к=0

2 (inverf (у0))2 V Нк (inverf (у0)).Z0 =

e ,

к=0 (к + 1)!

2 -х2^ Нк (х0).г0к+1

= у0-~TZe £ 0

к=0

Ур к=0 (к +1)!

Многочлены Эрмита могут быть вычислены на основе рекуррентных соотношений

Нп+1(х) = 2хНп (х) - 2пНп_1(х), где Но(х) =1Н1(х) = 2х.

Для уменьшения погрешности за счёт резкого возрастания порядков чисел ряда (1) перейдём к другому рекуррентному соотношению:

2

У к+1(х) = ~7—т(хП(х) _ У к _1(х)) (2)

к +1

где У0(х) = 1,^1(х) = 2х,...,Ук(х) =

Нк (х)

к!

В результате формула (1) может вычисляться на основе следующего рекуррентного соотношения:

к +1 z0

ук+1 = ук I—e °.ук(х0)., л . (3)

VP к +1

Для вычисления погрешности при оставлении в ряде (1) m членов следует подставить в невязку величину у0 = у + D0, где А0 - абсолютная погрешность приближения на заданном интервале. В результате получим [13]:

А(т) » 2 e~х2 Нm+1 (х0 ) ( eУ°2 ) m+1 Ат+1

0 4Р ' (т +1)^ 2 ^ 0 '

Начальные приближения целесообразно выбирать на основе сегментной (разновидность сплайн) аппроксимации. При использовании постоянных приближений на соответствующем подынтервале с учётом соответствующих норм погрешностей для х е [х1, хг+1] будет [13 ,16]:

у0,- =(erf (хг) +erf (хг+1 )V2;

у0, = erf (

либо у, = th( х).

При программной реализации алгоритмов на ЭВМ необходимая точность вычислений задается в виде e = 0,5.10~a, где степень а определяет количество верных значащих цифр результата. Тогда выбираем условие прекращения:

Ук - Ук-1

Ук

+

Ук+1 - Ук

Ук+1

=к-1

2 -X2 , . zi+1 /

Є 0gi СX0)^— Уг-

£ Є.

(4)

/Ж i +1/

Для расчета по формуле (3) задается начальное приближение (x0,у0) . В частности, можно задать (х0 = 0, у0 = 0). Для стабилизации времени счёта при любых значениях аргумента вводится система опорных точек. Диапазон изменения аргумента x, начиная с нуля, разбивается на интервалы с шагом h . В качестве опорной точки в каждом интервале берем его середину и вычисляем yi = erf (xi) ,

причем xt = h(i — 1) +1 h = h(i — 1). Из формул (3) и (4) видно, что количество шагов итераций зависит от h .

Для примера приведем две системы опорных точек (h = 0,5 и h = 0,25) и их соответствующие количества шагов итераций при программной реализации.

Таблица 1. Система опорных точек при h = 0,5

Таблица 2. Система опорных точек при h = 0,25

X0 Уо

0,25 0,276326390168236933

0,75 0,711155633653515132

1,25 0,922900128256458230

1,75 0,986671671219182444

2,25 0,998537283413318848

2,75 0,999899378077880363

3,25 0,999995697225363250

3,75 0,999999886272743430

4,25 0,999999998149425820

4,75 0,999999999981514951

X0 Уо X0 Уо

0,125 0,140316204801333839 2,625 0,999794624263858738

0,375 0,404116909434822325 2,875 0,999952145160256182

0,625 0,623240882188417924 3,125 0,999990103270201730

0,875 0,784075061059859613 3,375 0,999998184723399535

1,125 0,888388231701707819 3,625 0,999999704859807446

1,375 0,948170072782090378 3,875 0,999999957486055875

1,625 0,978443733239983643 4,125 0,999999994576599199

1,875 0,991990057670119949 4,375 0,999999999387516704

2,125 0,997345970640517631 4,625 0,999999999938783840

2,375 0,999217061782108854 4,875 0,999999999994586539

Рис. 1. График количества шагов итерации (линия соединения точек) при к = 0,5 и а = 7, 10, 15

Количество шагов итерации

14 13 12 11 10 9

3 7 6 5

4 3 2

1 2 3 4 5

Рис. 2. График количества шагов итерации (линия соединения точек) при к = 0,25 и а = 7, 10, 15 3. Вычисление тувг/ (х) и шуФ 0( х)

В работах [13, 16] приведены методы получения итерационных формул произвольного порядка сходимости на основе соответствующих функциональных уравнений и методов дифференцирования, что для получения итерационных формул высоких порядков представляется весьма трудоёмкой задачей. Итерационная формула второго порядка:

Ум = У! -^р(ег/(у.)- х).еуУ.

Итерационная формула третьего порядка:

У.+1 = У.-^р (ег/ (уг) - х).еУ2 +Р .у. (ег/ (уг) - х)2.е 2у2. (5)

Однако более просто получить эти формулы за счёт обращения ряда (1).

В основе обоснования такого типа обращений в соответствии с [17] рассмотрим решение уравнений рядами. Предположим, что функция у = / (х) представлена в неявном виде:

и ¥(х, у) в окрестности точки (х0, у0) разлагается в ряд по степеням х - х0 и у - у0 , причём постоянный член в нём равен нулю, а коэффициент при у - у0 отличен от нуля. Тогда и функция у = / (х), определяемся уравнением (6) в окрестности указанной точки, также разлагается в ряд по степеням х - х0 вблизи х = х0, т.е. если функция ¥(х,у) является аналитической в точке (х0,у0) , то и функция у = /(х) , определяемая уравнением (6), будет аналитической в точке х0.

Отсюда следует, что функция у = /(х) в некоторой окрестности точки х - х0 представляется

рядом:

где у0 - свободный член, то при a1 ф 0 в окрестности у = у0, из (7) величина х определяется как функция от у , разлагающаяся, в свою очередь, в ряд по степеням у — у0. Таким образом, если у является аналитической функцией от х в точке x0 , то в соответствующей точке y0 (при выполнении указанного условия) обратная функция будет аналитической).

Так как erf (x) является функцией аналитической, то можно обратить степенной ряд (1). Техника такого обращения хорошо освещена в работе [13]. В современных условиях для применения таких преобразований целесообразно использовать системы компьютерной алгебры [19 - 24].

В работах [23, 25] рассмотрен несколько иной подход, чем было указано выше. Так, вместо невязки z0 = x — x0 , удовлетворяющей функциональному уравнению F(x, у) = x — f (у) = 0 ,

используется невязка z0 = (x — f(у0))/f'^) , удовлетворяющая функциональному уравнению

F(x, у) = (x — f (у))/f'(у) = 0 , так что в конечном счёте итерационные формулы для обоих видов невязок совпадают с точностью перегруппировки членов. Однако для получения подобных формул уже используются прямые методы разложения по невязкам. В этом случае для inverf (x) невязка имеет

Очевидно, такой тип невязок наиболее подходит при вычислении интегральных функций. Используя первый вид невязок, т.е. z0 = x — erf (у0) на основе обращения ряда невязок (1) (таблица) и подставляя значения полиномов Эрмита и невязки получим

F С x, у) = О

(6)

у - Уо = a1 С x - Xo) + a2Сx - Xo)2 +... + an ^ - Xo)n +...,

(7)

Р 2

у = іпуєг/(х) = ітег/(ег/(у0) - г0) = у0 - -у-еу° .г0 +

2у0(^Р еу2°.^0)2 2(4у0 + 1)^~Р еу2°.^0)3

+--------2----------------------------2-----------+... .

2! 3!

\/ х 2

Для невязки 20 =—^~ (х - ег/ (у0)).еу° получим [23]

2 /1 4 2 \ 3

у = у0 + 20 + у0.^0 + (3 + 3.у0).^0 +... .

Погрешности при оставлении двух, трёх и т.д. членов составляют [13, 23]

А2 »-у.А2 -(3у2 -|).А3;

Аз »(3у2 + 3).А3 + -1(4у2 + 5).у.А4;

А4 »-6(12у2 + 7).у.А4 + (24у4 -358у2 - 15-).А5.

6 5 5 15

Отметим, что погрешности для первого и второго подходов выбора невязок совпадают.

Для сравнения полученных результатов с известными приведём результаты алгоритма

Н.Е. Рейв [13] для вычисления у = туег/с( х) по формулам

у0 =-(1п^л/р(- 1п х)1/2));

у п+1 =(- 1п ¥ (х уп ))1/2;

¥ (у, х) = ^ 3^

1 + 1 + 1 + 1 +...

где г = 1/ у.

Необходимое количество членов непрерывной дроби ! и итераций к при точности 12 десятичных знаков:

х 10-6 10 -4 10 -2 0,05 0,1 0,2 0,3

і 29 33 51 74 98 150 220

к 7 8 11 14 26 21 27

Отметим, что для значительно большей разрядности, достигающей 30 десятичных разрядов, при счёте по предложенным итерационным формулам количество итераций колеблется в диапазоне 4

- 7 для всего интервала изменения аргумента х є [0,5], что ведёт к значительно меньшему времени счёта по сравнению с вышеприведенным алгоритмом.

При программной реализации итерационной формулы (5) сходимость вычислений тесно связана с начальным приближением. Авторы предлагают следующее: для маленьких значений аргумента, примерно x < 0,2 , в качестве начального приближения следует выбирать

логарифмическую формулу у0 = -д/- 1п(1 - x2) , а для остальных значений аргумента, воспользуемся линейной сегментной аппроксимацией: у0 = C^ -1) + D , где коэффициенты С и й определяются в зависимости от интервала изменения аргумента x.

Таблица 3. Коэффициенты линейной сегментной аппроксимации функции ітег^х)

Интервал изменения аргумента х С й

[0;1 -10-2) 1,8398 1,5559

[1 -10-2;1 -10 -4) 93,902 2,5665

[1 -10-4;1 -10 -6) 7149,4 3,3124

-1 о і -1 о ! оо 599310 3,9270

[1 -10-8;1 -10-10) 5,2592.107 4,4620

[1 -10-10;1 -10-12) 4,7147.109 4,9399

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

распределения и соответствующих квантилей может быть использован не только для вычисления

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

Помимо этого, описанный в статье подход для вычисления интеграла вероятностей и

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

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

которому с заданной вероятностью Р соответствует условие х < t , т.е. выполняется равенство

г

| / (Х)^Х = Р, где / ( х) - функция плотности распределения.

В нашем случае квантиль t адекватен функциям inverf (х) либо invF0(х).

4. Проверка согласованности результатов вычислений

Для проверки согласованности вычисляем величины:

А пр = erf (inverf (х)) - х и А обр = inverf (erf (х)) - х при различных значениях аргумента х и

точности вычислений a . Напомним, что £ = 0,5.10-a. Результаты такой проверки приведены на рис. 3 и 4

Testing Module

X

erf(inverf(X)) -X

0.001 0.000000000000000 0.002 0.000000000000000 0.003 0.000000000000000

0.004 0.000000000000000

0.005 0.000000000000000

0.006 0.000000000000000 0.007 0.000000000000000

0.008 0.000000000000000 0.00Э 0.000000000000000

0.01 0.000000000000000 0.011 0.000000000000000

0.012 0.000000000000000

а) при a = 9 и х » 0 б) при a = 9 и х » 1

Рис. 3. Результаты вычислений величины А п

|| _ Testing Module

X erf(inverf|X)) - X

0.988 0.000000000001357

0.989 0.000000000002102

0.99 -0.000000000000898

0.991 -0.000000000000807

0.992 0.000000000000621

0.993 -0.000000000000607

0.994 -0.000000000001012

0.995 -0.000000000001733

0.996 0.000000000004420

0.997 0.000000000000507

0.998 -0.000000000003556

0.999 0.000000000004338

Testing Module

X erf(inverf|X)) - X

0.988 0.000000000000000

0.989 0.000000000000000

0.99 0.000000000000001

0.991 0.000000000000001

0.992 -0.000000000000002

0.993 0.000000000000000

0.994 0.000000000000010

0.995 -0.000000000000006

0.996 -0.000000000000002

0.997 0.000000000000003

0.998 -0.000000000000001

0.999 0.000000000000002

в) при a = 12 и х » 1 : erf (inverf (х)) - х

I Testing Module

X inverf(erf(X)) - X

ояП 0.000000000000

0.02 0.000000000000

0.03 0.000000000000

0.04 0.000000000000

0.05 0.000000000000

0.06 0.000000000000

0.07 0.000000000000

0.08 0.000000000000

0.09 0.000000000000

0.1 0.000000000000

0.11 0.000000000000

0.12 0.000000000000

I Testing Module

X inverf(erf(X)) - X

3.89 0.000000001751

3.9 0.000000001302

3.91 -0.000000000738

3.92 0.000000005500

3.93 0.000000006855

3.94 0.000000008909

3.95 0.000000004610

3.96 0.000000003627

3.97 -0.000000000979

3.98 -0.000000003784

3.99 0.000000003438

14 0.000000006770

Testing Module

X inverf(erf(X)) - X

3.89 0.000000000002

3.9 0.000000000000

3.91 0.000000000000

3.92 -0.000000000002

3.93 -0.000000000006

3.94 -0.000000000014

3.95 -0.000000000033

3.96 0.000000000008

3.97 0.000000000013

3.98 0.000000000028

3.99 0.000000000003

14 0.000000000002

а) при a = 12 и х » 0 б) при a = 12 и х » 1

Рис. 4. Результаты вычислений величину А обр

в) при a = 15 и х » 1 : inverf (erf (х)) - х

5. Выводы

В статье приведены алгоритмы вычисления функции интеграла вероятностей для нормального распределения и квантилей на основе разложения этих функций в ряд невязок и приведены результаты их программной реализации.

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

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

СПИСОК ЛИТЕРАТУРЫ

1. Справочник по специальным функциям с формулами, графиками и математическими таблицами I Под ред. М. Абрамовица и И. Стиган: Пер. с англ. - М.: Наука, 1979. - S32 с.

2. Градштейн И.С. , Рыжик И.М. Таблицы интегралов, сумм, рядов и произведений. - М.: Физматгиз, 1963. - 1100с.

3. Корн Г., Корн Т. Справочник по математике для научных работников и инженеров. - М.: Физматгиз, 197S. - 832с.

4. Бронштейн И.Н., Семендяев К.А. Справочник по математике для инженеров и учащихся втузов. - М.: Наука, 19S1. - 72G с.

б. Янке Е., Эмде Ф., Лёш Ф. Специальные функции. Формулы, графики, таблицы. - М.: Наука, 1964. - 344с.

6. Прудников А.П., Брычков Ю.А., Марычев О.И. Интегралы и ряды. Специальные функции. - М.: Наука, 19S3. -7б2 с.

7. Таблицы нормального интеграла вероятностей, нормальной плотности и её нормированных производных I Под ред. чл. - корр. АН СССР Н.В. Смирнова. - М.: Издательство АН СССР, 196G. - 136 с.

5. Большев Л.Н., Смирнов Н.В. Таблицы математической статистики. - М.: Наука, 19S3. - 416 с.

9. Хаджи П.И. Функция вероятности. Интегралы, ряды и некоторые обобщения. - Кишинев: Издательство АН Молдавской ССР, 1971. - 39S с.

1G. Fike C.T. Computer evaluation of mathematical function. - New Jersey: Prentice - Hall, 196S. - 22S p.

11. Computer approximations. - Ed. by J.F. Hart. New York: Wiley, 196S. - 344 p.

12. Люк Ю. Специальные математические функции и их аппроксимации. - М.: Мир. - 19SG. - 6GS с.

13. Попов Б.А. , Теслер Г.С. Вычисление функций на ЭВМ: Справочник. - Киев: Наукова думка, 19S4. - 6GG с.

14. Philip J.R. The function invert^. “ Austral. J. Phys.” 13. - 196G. - P. 13 - 2G.

16. Strecok A.J. On the calculation of the inverse of the error function. “ Math. Comput.” 22. - 196S. - P. 144 - 158.

16. Попов Б.А., Теслер Г.С. Приближение функций для технических приближений. - Киев: Наукова думка, 19SG. -366 с.

17. Фихтенгольц Г.М. Курс дифференциального и интегрального исчисления. - М.: Физматлит, 1962. - SG7 с. (С. б01).

18. Фильчаков П.Ф. Численные и графические методы прикладной математики: Справочник. - Киев: Наукова думка, 197G. - SGG с.

19. Heck A. Introduction to Maple. - N.Y. : Springer - Verlag, 1993. - 498 p.

2G. Maeder R.E. Programming in Mathematica. - Redwood City, California: Addison - Weslagerg, 199G. - 267 p.

21. Hearn A.C. REDUCE User Manual. - Santa Monica: The Rand Corporation, 1996. - 480 p.

22. Misniewski M., Zero A. MathCAD Plus 6.G. - Warszawa: EXIT, 1998. - 322 p.

23. Попов Б.О. Розв'язування математичних задач у системах комп'ютерної алгебри. - MAPLEV. - Київ: V i P, 2GG1.

- 312 с.

24. Морозов А.А., Клименко В.П., Фишман Ю.С. АНАЛИТИК - 2GGG II Математические машины и системы. - 2GG1. № 1. - C. 66 - 99.

26. Попов Б.О., Лушник О.І. Побудова оптимального алгоритму для обчислення функції, оберненої до інтегралу імовірності II Волинський мат. вісник. - 1999. - Вип. 6. - C. 1G7 - 112.

26. Митропольский А.К. Интеграл вероятностей. - Л.: Изд-во Ленинградского университета, 1972. - 86 с.