Если ввести формулу =NORM.DIST(115,100,15,TRUE) в ячейку, Excel без лишних сложностей вернет значение 0.8413447. Этот вызов похож на простой поиск по таблице готовых значений. Однако это не так. За этим числом скрывается интегральная функция нормального распределения, не имеющая аналитического выражения в элементарных функциях, а функции CHISQ.INV.RT и BETA.DIST базируются на вычислении специальных математических функций, которые библиотека должна рассчитывать с высокой точностью, а не подбирать приближенно. Табличный компонент, заявляющий о совместимости с Excel, должен воспроизводить эти значения до последнего знака, что требует точной реализации численных методов, а не просто дублирования имен функций.
HotXLS поддерживает более пятидесяти таких статистических функций, и работа по обеспечению их точности полностью скрыта от глаз пользователя формул. В этом обзоре рассматриваются принципы расчетов движка: общее математическое ядро специальных функций, ветвление вычислений для сохранения их стабильности, а также одна ошибка расчета обратного нормального распределения в крайних значениях диапазона, которая долго оставалась незамеченной, так как обычные тесты ее не затрагивали.
Один вызов на листе и пятьдесят распределений под капотом
Реализованные функции покрывают все основные типы распределений, применяемые при статистическом анализе. Сюда входят нормальное распределение (NORM.DIST и NORM.S.DIST с их обратными функциями), семейство гамма- и хи-квадрат распределений (GAMMA.DIST, CHISQ.DIST, CHISQ.DIST.RT, CHISQ.INV.RT), бета-распределение (BETA.DIST и BETA.INV), распределения выборок (T.DIST, T.DIST.2T, F.DIST и F.INV), дискретные распределения (BINOM.DIST и POISSON.DIST), а также вспомогательные функции доверительных интервалов (например, CONFIDENCE.T и CONFIDENCE.NORM). С точки зрения разработчика любая из них представляет собой простую формулу. Вы задаете входные данные в ячейках, запускаете вычисление книги и считываете готовый результат.
var
wb: IXLSWorkbook;
sh: IXLSWorksheet;
begin
wb := TXLSWorkbook.Create;
sh := wb.Sheets.Add;
sh.Range['A1', 'A1'].Value := 115; // observation
sh.Range['A2', 'A2'].Value := 100; // mean
sh.Range['A3', 'A3'].Value := 15; // standard deviation
// The XLS formula parser uses ';' as the argument separator.
Writeln(wb.Calculate('=NORM.DIST(A1;A2;A3;TRUE())')); // 0.8413447
Writeln(wb.Calculate('=CHISQ.INV.RT(0.05;10)')); // 18.3070381
Writeln(wb.Calculate('=BETA.DIST(0.5;2;3;TRUE())')); // 0.6875
end;Метод Calculate книги компилирует и вычисляет формулу «на лету» в контексте текущего листа, возвращая значение типа Variant. Важно помнить об одной детали: парсер формул метода Calculate использует точку с запятой в качестве разделителя аргументов, то есть синтаксис имеет вид =SUM(A1;B1), а не =SUM(A1,B1). При этом формулы, сохраненные непосредственно в ячейках листа, используют стандартную для Excel запятую. Этот же вычислитель обрабатывает и все статистические функции, поэтому при успешной работе любого вызова в Calculate все остальные функции вычисляются по тому же пути.
Две базовые функции как основа всех расчетов
Большинство функций интегрального распределения из этого набора не рассчитываются прямым суммированием или интегрированием их формул. Вместо этого вычисления строятся на основе двух специальных функций: регуляризованной неполной гамма-функции P(a, x) и регуляризованной неполной бета-функции Ix(a, b). Внутренняя структура вызовов довольно проста. Функция распределения хи-квадрат представляет собой функцию гамма-распределения с параметром формы df/2 и масштабом 2. Функция гамма-распределения рассчитывается напрямую через P(a, x). Интегральные функции распределений Стьюдента (t), Фишера (F) и биномиального распределения выражаются через значения регуляризованной неполной бета-функции при соответствующих аргументах. Функция распределения Пуассона вычисляется через верхнюю неполную гамма-функцию Q. Качественная реализация гамма- и бета-функций автоматически гарантирует точность вычисления десятка других распределений.
Понятие «регуляризованная» играет здесь ключевую роль. Обычная неполная гамма-функция растет со скоростью факториала, а интеграл бета-функции может приводить к потере значимости или переполнению задолго до получения конечного ответа. Регуляризованные формы делятся на полную гамма- или бета-функцию, благодаря чему их значения всегда лежат в интервале от нуля до единицы, что точно соответствует диапазону вероятностей. Эта нормализация позволяет одному и тому же коду вычислять распределение хи-квадрат как для двух, так и для двухсот степеней свободы без риска выхода промежуточных значений за границы диапазона double. Это также объясняет, почему функцию распределения нельзя вычислять простым сложением длинного хвоста плотностей: каждый член вносит погрешность округления, эти погрешности накапливаются в процессе суммирования, в то время как регуляризованная специальная функция заменяет сумму вычислением быстро сходящегося ряда или непрерывной дроби.
Степенной ряд ниже диагонали и непрерывная дробь выше нее
Перед началом вычислений неполной гамма-функции движок принимает одно решение: сравнивает x со значением a + 1. Этот порог выбран не случайно. Степенной ряд для P(a, x) быстро сходится, если значение x мало относительно a, и сходится крайне медленно или неэффективно при больших x. Непрерывная дробь демонстрирует противоположные свойства. Движок использует степенной ряд при x меньше a + 1 и непрерывную дробь Ленца при x больше или равном a + 1, выполняя в каждой ветви наиболее эффективный для данных условий расчет.
Вычисление непрерывной дроби требует одной проверки защиты. Метод Ленца работает с накоплением числителя и знаменателя с инвертированием знаменателя на каждом шаге, и если любое из этих значений приближается к нулю, операция деления вызывает сбой. Решением является минимальный порог: если промежуточный член становится по модулю меньше примерно 1e-30, его значение принудительно приравнивается к 1e-30, что гарантирует конечность рекуррентного соотношения без искажения сходящегося результата. По этой же причине аналогичное ограничение применяется и в непрерывной дроби для неполной бета-функции. Эта небольшая константа выполняет важную стабилизирующую роль, предотвращая деление на ноль.
Верхний хвост Q(a, x) рассчитывается просто как 1 минус P(a, x), и именно так вычисляется интегральная функция распределения Пуассона: вероятность регистрации не более k событий при среднем λ равна Q(k + 1, λ). Использование верхней неполной гамма-функции вместо сложения k + 1 членов распределения Пуассона является более рациональным выбором в пользу расчета одного сходящегося выражения.
Вычисление дискретных вероятностей без переполнения факториала
Дискретные распределения несут в себе иную опасность. Расчет биномиальной вероятности включает вычисление биномиального коэффициента, и, например, число сочетаний из 52 по 26 представляет собой гигантское число. Прямой расчет приведет к переполнению типа double в числителе до выполнения деления, возвращающего значение к диапазону вероятностей. Движок никогда не выполняет вычисления напрямую. Он рассчитывает факториалы в логарифмическом пространстве с помощью логарифма гамма-функции, складывает и вычитает полученные логарифмы, добавляет логарифмы вероятностей успеха и неудачи, и только в самом конце выполняет возведение в степень (экспоненту).
// Binomial probability mass, evaluated entirely in log space.
// LnGammaF(n+1) is ln(n!)
result := Exp(LnGammaF(nt + 1) - LnGammaF(kk + 1) - LnGammaF(nt - kk + 1)
+ kk * Ln(pp) + (nt - kk) * Ln(1 - pp));Сама логарифмическая гамма-функция представляет собой аппроксимацию Ланцоша, точную во всем положительном диапазоне и не требующую больших вычислительных ресурсов. Поскольку все крупные промежуточные значения хранятся в виде логарифмов до финального вызова Exp, наибольшее число, возникающее при вычислениях, не превышает единицу (сама вероятность). Функция плотности распределения Пуассона вычисляется по тому же алгоритму, где один логарифм гамма-функции заменяет факториал в знаменателе. Граничные случаи (когда вероятность p равна точно нулю или единице) обрабатываются отдельно, чтобы код не вызывал функцию Ln(0). HotXLS возвращает 0.2460938 для вызова BINOM.DIST(5,10,0.5,FALSE) and 0.6766764 для интегральной функции POISSON.DIST(2,2,TRUE), что полностью совпадает с результатами расчетов в Excel.
Вычисление обратных функций методом локализации корня
Обратная функция распределения решает обратную задачу: по заданной вероятности найти значение, интегральная функция распределения которого равна этой вероятности. Только одна обратная функция из этого набора имеет простую прямую формулу. Функция NORM.S.INV (обратное стандартное нормальное распределение) использует рациональную аппроксимацию Аклама - соотношение двух полиномов, обеспечивающее точность на уровне типа double во всем диапазоне и разделенное на центральную область и два хвоста распределения. Это вычисление в закрытой форме, не требующее итерационного подбора.
Другие обратные функции не имеют аналитических формул, поэтому движок вычисляет их численно. Он ограничивает искомое значение нижним и верхним пределами из области определения распределения, а затем выполняет деление пополам: вычисляет прямую функцию распределения в средней точке, сдвигает границу интервала для удержания целевой вероятности и повторяет процесс до достижения нужной точности. Для обратных функций гамма-распределения и распределения хи-квадрат границы начинаются от нуля до примерной верхней оценки, которая удваивается, если вероятность не попадает в интервал. Обратная функция распределения Стьюдента (t) использует симметричные расширяющиеся границы, а обратное распределение Фишера (F) выполняет деление пополам на неотрицательном интервале. Это требует всего нескольких десятков вычислений прямой функции распределения, что незаметно при расчетах, но гарантирует, что обратная функция имеет ту же точность, что и прямая. Именно поэтому круг операций вроде CHISQ.DIST(CHISQ.INV(0.7,5),5,TRUE) с высокой точностью возвращает 0.7.
Десятичный логарифм, скрывавшийся в хвосте распределения
В истории разработки был один примечательный баг, способный долго оставаться незамеченным. Алгоритм Аклама для вычисления обратного нормального распределения содержит три ветви вычислений. Широкая центральная ветвь, применяемая при значениях вероятности от 0.025 до 0.975, вычисляет отношение полиномов без использования логарифмов. Две боковые ветви (хвосты распределения) для очень малых или очень больших вероятностей сначала вычисляют логарифм входного значения, так как поведение хвоста аппроксимируется квадратным корнем из отрицательного натурального логарифма вероятности p.
В ранней версии кода в ветви хвостов распределения ошибочно использовался десятичный логарифм вместо натурального. Разница между ними составляет константный коэффициент примерно 2.30, поэтому результаты расчетов в хвостах стабильно расходились на заметную величину. Тем не менее, при обычных проверках функция работала корректно, так как большинство тестов выполнялись в центральном диапазоне. Вызовы NORM.S.INV(0.5) возвращают 0, NORM.S.INV(0.975) возвращают табличное значение 1.959964, и оба эти расчета идут по центральной ветви без логарифмов. Ошибка проявлялась лишь тогда, когда вероятность переходила в хвост распределения (например, при NORM.S.INV(0.001), которое должно возвращать -3.0902323). Любые зависимые от обратного нормального распределения функции в хвостах унаследовали это смещение. Вывод прост и очевиден: функции с ветвлением алгоритма требуют тестирования каждой ветви по отдельности, так как корректная работа основной логики маскирует сбои в редких сценариях. Исправление свелось к замене функции десятичного логарифма на натуральный, после чего крайние значения стали точно совпадать с результатами Excel.
Знак аргумента x определяет направление хвоста t-распределения
Интегральная функция распределения Стьюдента (t) содержит один нюанс, который легко перепутать. Ее значение рассчитывается на основе регуляризованной неполной бета-функции со значением аргумента df / (df + x²), но результат бета-функции отражает вероятность в хвосте за пределами модуля x, а не интегральную вероятность до точки x. Из-за симметричной формы t-распределения алгоритм пересчета зависит от знака аргумента x.
// Student t CDF. ib is the regularized incomplete beta at df/(df+x*x),
// which measures the symmetric tail. The cumulative value depends on
// the sign of x; returning ib unconverted gives the wrong tail.
ib := BetaIF(df / 2, 0.5, df / (df + x * x));
if x > 0 then
result := 1 - 0.5 * ib // above the mean: one minus half the tail
else if x < 0 then
result := 0.5 * ib // below the mean: half the tail
else
result := 0.5; // exactly at the meanПри x больше нуля интегральная вероятность равна единице минус половина симметричного хвоста; при x меньше нуля она равна половине этого хвоста; в нуле она составляет ровно 0.5. Если просто возвращать значение бета-функции напрямую, результат будет соответствовать противоположной стороне распределения с погрешностью на величину всей кривой для любого ненулевого x. Варианты для правого хвоста и двустороннего распределения строятся на базе этого же алгоритма, поэтому вызов T.DIST.2T(1,1) возвращает 0.5, а T.DIST(1,1,TRUE) возвращает 0.75, а расчет обратной функции T.INV использует метод деления пополам для точного совпадения значений.
Все эти сложные расчеты скрыты от пользователя, что и является главной целью. Вы просто вводите формулу и получаете результат, совпадающий с Excel. Если вы хотите дополнить движок собственной бизнес-логикой, правила регистрации функций описаны в нашем обзоре движка формул и пользовательских функций, а механизм работы со ссылками между листами и именованными диапазонами представлен в статье об именованных диапазонах и межстраничных формулах. Все эти возможности входят в состав библиотеки HotXLS spreadsheet component для Delphi и C++Builder наряду с API для чтения, записи, построения диаграмм и форматирования таблиц.