Как Excel считает функцию нормального распределения?

[faruk]

Как Excel считает неберущийся интеграл для функции нормального распределения?

$\Phi(x) = \frac{1}{\sqrt{2\pi}} \int\limits_{-\infty}^x e^{-\frac{t^2}{2}} \, dt$

[Sonic86]

В принципе он может разбивать $\int\limits_{-\infty}^x=\int\limits_{-\infty}^0+\int\limits_0^x$, первый интеграл ясно чему равен, а второй интеграл - через ряд Маклорена (он все равно быстро сходится). Это чисто теоретически. Как на самом деле - не знаю.

[AlexValk]

Для вычисления специальных функций численно (в том числе в комплексной области) обычно используют специальные быстрые алгоритмы. Подробнее посмотрите, например книгу
"Люк Ю. Специальные математические функции и их аппроксимации."
Общая идея (неплохо работающая для функции $\Phi(x)$):
1. Для больших по модулю $x$ использовать асимптотические разложения.
2. Для небольших по модулю $x$ использовать ряды Тейлора и (то, что указал Sonic86) или некоторые другие разложения (например, апроксиманты Паде).
3. В эти формулы вводятся небольшие коррекции для лушей склейки "склейки" этих формул в области промежуточных $x$.
4. Бывает удобно предварительное преобразование переменной.
Явные формулы для $\Phi(x)$) (вплоть до точностей $\sim 10^{-20}$) можно посмотреть в книге Люка, а при меньших запросах к точности, например в книге Абрамовиц, Стиган "Справочник по специальным функциям" (и во многих других, посвященных численным расчетам). (Обычно в книгах даются формулы для ${\rm erf}(x)$ - функции ошибок, которая тривиально связана с $\Phi(x)$)

Sonic86 в сообщении #568079 писал(а):

Как на самом деле - не знаю.

Какой именно алгоритм вычисления $\Phi(x)$ используют в Excel программисты Microsoft (и меняется ли он от версии к версии) я тоже не знаю, хотя предполагаю - стандартный (из книги Люка). А вот, что точно можно сказать - алгоритм, описанный здесь:

Sonic86 в сообщении #568079 писал(а):

В принципе он может разбивать $\int\limits_{-\infty}^x=\int\limits_{-\infty}^0+\int\limits_0^x$, первый интеграл ясно чему равен, а второй интеграл - через ряд Маклорена (он все равно быстро сходится).

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

$$\Phi(x) =\frac12+\frac{x}{\sqrt{2\pi}}\sum_{k=0}^\infty\frac{(-1)^k}{2^k k!(2k+1)}x^{2k}.$$

Этот ряд действительно сходится для всех $x\in \mathbb{R}$, но хорош он для численного использования только при небольших $|x|$. Если его попытаться использовать для вычисления, например, хотя бы величины $\Phi(20)$, то возникнут практически непреодолимые численные трудности. Дело в том, что поначалу члены ряда быстро растут по модулю и лишь после этого начинают убывать, стремясь в итоге к нулю. Например, 200-й член ряда для $\Phi(20)$ превышает $10^{84}$ и лишь 536-й член становится $<1$. В такой ситуации необходимо будет проводить вычисления со 100 значащими цифрами, чтобы после сокращения колоссальных положительных и отрицательных слагаемых получить правильный ответ (хотя бы из таблиц, хорошо известно, что с очень высокой точностью $\Phi(20)=1$).
Похожие проблемы есть не только для специальных функций, но и для элементарных - например, если бы мы попытались вычислить $\sin(x)$ по его тейлоровскому ряду при достаточно больших $x$. На практике в этом случае нормальный численный алгоритм сначала, используя периодичность $\sin$ "сбросит" $x$ в интервал $[-\pi;\pi]$ (а еще лучше - используя формулы приведения - в интервал $[-\pi/2;\pi/2]$), и уже только потом применит Тейлоровский ряд (или еще что-то подобное). Так, что без использования "аналитики" во многих случаях не удастся создать эффективный численный алгоритм. На заре компьютерной эры эти вопросы активно разрабатывались и было предложено немало достаточно эффективных численных алгоритмов вычисления элементарных и специальных функций, которые теперь реализованы в инженерных калькуляторах, стандартных подпрограммах языков программирования и математических программах. По этой причине, кстати, к результатам численных расчетов значений спецфункций в математических программах (даже таких "монстров" таких как Mathematica или Maple), особенно в комплексной области, следует относится с осторожностью - всегда не лишне их проверять на "здравый смысл": нет 100% гарантии, что используемый данной программой численный алгоритм правильно обрабатывает все ситуации.

[Sonic86]

AlexValk, спасибо, я про асимптотические разложения забыл

[PAV]

Учитывая наличие хороших содержательных ответов, тема перенесена из "свободного полета" в тематический раздел "численные методы"

[AlexValk]

Добавлю, что в настоящее время многие классические численные алгоритмы (в том числе алгоритмы вычисления элементарных и специальных функций) требуют пересмотра (и реально пересматриваются) для эффективного использования многоядерности и большой оперативной (и кэш) памяти современных компьютеров. Поэтому нет гарантии, что от версии к версии Excel алгоритмы не меняются.
Хотя быстрота численных алгоритмов вычисления спецфункций, вероятно, совсем не первая вещь, которую стремятся оптимизировать в Excel - большинству пользователей Excel более важны общая эффективность использования памяти, компактность форматов данных для быстрой записи-чтения, и т.д.. За численной эффективностью больше следят производители математических программ (таких как Mathematica, Maple, Matlab и т.д.) и библиотек подпрограмм численных расчётов (таких как NAG - Numerical Algorithms Group, Visual Numerics IMSL и т.п.). Впрочем недавно Microsoft объявила библиотеку Cloud Numerics для платформы .NET - чтобы эффективно выполнять многопроцесорные облачные вычисления.
Отдельная область современных численных алгоритмов связана с использованием GPU (вместо CPU) для численных вычислений, производительность которых в этих задачах часто в десятки раз быстрее чем у лучших CPU.

[ewert]

AlexValk в сообщении #568187 писал(а):

уже только потом применит Тейлоровский ряд (или еще что-то подобное).

Тейлоровский не стоит -- у него сходимость всё-таки не ахти (если аргумент не мал). Лучше применять чебышёвское приближение.

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