Введіть =NORM.DIST(115,100,15,TRUE) у клітинку, і Excel просто поверне 0.8413447. Цей виклик виглядає як звичайний пошук за таблицею. Але це не так. За цим одним числом стоїть інтегральна функція нормального розподілу, яка не має закритої аналітичної форми, а за CHISQ.INV.RT та BETA.DIST стоять спеціальні математичні функції, які бібліотека повинна точно обчислювати, а не грубо оцінювати вручну. Компонент електронних таблиць, який претендує на повну сумісність з Excel, повинен відтворювати ці значення до останнього знака, який показує Excel, що означає відтворення чисельних методів, а не просто назв функцій.
HotXLS реалізує понад п'ятдесят таких статистичних функцій, і обчислювальна робота, яка робить їх точними, майже повністю прихована від рядка формул. Це огляд того, як двигун розраховує ці функції: спільне ядро спеціальних функцій, вибір гілок обчислень для збереження стабільності безпеки математичних операцій та одна помилка оберненого нормального розподілу, яка довгий час ховалася в 'хвості' розподілу, оскільки звичайні випадки ніколи не зачіпали цю проблемну ділянку коду.
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, інші слідуватимуть цим же шляхом.
Дві функції, на яких побудовано все інше
Більшість інтегральних функцій розподілу (CDF) у цьому наборі обчислюються не шляхом підсумовування або інтегрування їхніх власних визначень. Вони обчислюються на основі двох спеціальних функцій: регуляризованої нижньої неповної гамма-функції, яка записується як P(a, x), та регуляризованої неповної бета-функції, яка записується як Ix(a, b). Внутрішньо це помічники, на які спирається двигун, і ланцюжок зв'язків тут короткий. CDF розподілу хі-квадрат — це CDF гамма-розподілу з формою df/2 та масштабом 2. CDF гамма-розподілу — це безпосередньо P(a, x). Інтегральні функції розподілів Стьюдента (t), Фішера (F) та біноміального розподілу є значеннями регуляризованої неповної бета-функції при відповідних аргументах. CDF розподілу Пуассона — це верхня неповна гамма-функція Q. Якісна реалізація гамма- та бета-функцій дозволяє автоматично отримати високу точність для десятка розподілів.
Слово 'регуляризована' тут є ключовим. Звичайна неповна гамма-функція зростає подібно до факторіалу, а неповний інтеграл бета-функції може призвести до збою через втрату значущості або переповнення задовго до отримання результату. Регуляризовані форми діляться на повну гамма- або бета-функцію, тому їхні значення лежать виключно в діапазоні від нуля до одиниці, що точно відповідає діапазону ймовірностей. Ця нормалізація дозволяє одній і тій самій процедурі обслуговувати розподіл хі-квадрат як з двома ступенями свободи, так і з двомастами, без ризику виходу проміжних значень за межі типу double. Це також пояснює, чому не варто обчислювати CDF шляхом складання якогось довгого хвоста значень щільності: кожен елемент додає свою помилку округлення, ці помилки накопичуються під час підсумовування ряду, а регуляризована спеціальна функція дозволяє уникнути цього підсумовування шляхом обчислення ряду, що швидко збігається, або ланцюгового дробу.
Ряди під діагоналлю, ланцюговий дріб над нею
Процедура обчислення неповної гамма-функції приймає одне рішення перед початком будь-яких обчислень: вона порівнює 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 ще до виконання ділення, яке мало б повернути його до реального значення ймовірності. Двигун ніколи не обчислює його напряму. Він розраховує факторіали в логарифмічному просторі за допомогою логарифмічної гамма-функції (log-gamma), додає та віднімає логарифми, додає логарифми ймовірностей успіху та невдачі та виконує піднесення до експоненти лише один раз наприкінці.
// Binomial probability mass, evaluated entirely in log space.
// LnGammaF(n+1) is ln(n!); the three log-factorials form ln(C(n,k)),
// and the whole exponent is built before a single Exp call.
// ln P(X=k) = ln(n!) - ln(k!) - ln((n-k)!) + k*ln(p) + (n-k)*ln(1-p)
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) та 0.6766764 для інтегральної POISSON.DIST(2,2,TRUE), що точно відповідає значенням Excel.
Обернені функції шляхом взяття у дужки прямої кривої
Лише одна обернена функція в цьому наборі має швидку пряму формулу. NORM.S.INV, обернений стандартний нормальний розподіл, використовує раціональне наближення Аклама — пару поліноміальних відношень, точність яких наближається до точності типу double у всьому діапазоні, розділеному на центральну область та два хвости. Це аналітичне обчислення, яке не потребує ітерацій.
Інші обернені функції не мають подібної формули, тому двигун обчислює їх чисельно. Він обмежує шукану відповідь нижньою та верхньою межами, вибраними з області визначення розподілу, а потім застосовує метод дихотомії (ділення навпіл): обчислює пряму CDF у середній точці, зміщує ту межу, яка залишає цільову ймовірність всередині інтервалу, і повторює процес до звуження інтервалу. Для обернених функцій гамма-розподілу та хі-квадрат цей інтервал починається від нуля та щедрого верхнього оціночного значення, побудованого на основі форми та масштабу розподілу, подвоюючи верхню межу, якщо ймовірність ще не потрапила в інтервал. Обернена функція розподілу Стьюдента (t) обмежує симетричні межі, які розширюються назовні; обернена функція розподілу Фішера (F) ділить навпіл невід'ємний інтервал. Вартість цього — кілька десятків обчислень CDF на виклик, що є непомітним для швидкості роботи електронних таблиць, а перевага полягає в тому, що кожна обернена функція має точно таку ж точність, як і пряма функція, яку вона обертає. Ось чому зворотне перетворення на кшталт CHISQ.DIST(CHISQ.INV(0.7,5),5,TRUE) повертає значення 0.7 з високою точністю.
Десятковий логарифм, який сховався у хвості
У ранній версії хвостової гілки замість натурального логарифма використовувався десятковий логарифм. Ці два логарифми відрізняються на постійний множник приблизно в 2.30 разу, тому результати на хвостах розподілу містили системну і значну помилку. Проте при звичайних перевірках функція виглядала правильною, оскільки більшість перевірок стосуються саме середньої зони. Значення NORM.S.INV(0.5) дорівнює нулю, 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 менше нуля вона дорівнює половині цього хвоста; у нулі вона становить рівно половину. Якщо повернути значення бета-функції безпосередньо, ви отримаєте ймовірність не з того боку розподілу, яка відрізнятиметься на всю площу під кривою для будь-якого ненульового x. Варіанти з правим хвостом та двома хвостами будуються на цій самій гілці, саме тому T.DIST.2T(1,1) повертає 0.5, а T.DIST(1,1,TRUE) — 0.75, а обернена функція T.INV виконує пошук за цією скоригованою CDF, що забезпечує правильність зворотного перетворення.
Ніщо з цього не видно з клітинки таблиці, і саме такий результат і планувався. Ви пишете формулу та отримуєте число, яке збігається з Excel. Якщо ви розширюєте двигун власною логікою, принципи реєстрації функцій описані в нашій інструкції з двигуна формул та користувацьких функцій, а спосіб взаємодії формул між різними аркушами та іменованими діапазонами — в статті про визначені імена та формули між аркушами. Усе це постачається у складі компонента електронних таблиць HotXLS для Delphi та C++Builder разом з API читання, запису, діаграм та форматування, що розглядаються в інших публікаціях цього блогу.