Въведете =NORM.DIST(115,100,15,TRUE) в някоя клетка и Excel ще върне 0.8413447 без никакви формалности. Това извикване изглежда като обикновено търсене на стойност, но не е. Зад това единично число стои кумулативното нормално разпределение — интеграл, който няма аналитичен вид (closed form), а зад CHISQ.INV.RT и BETA.DIST стоят специални функции, които една сериозна библиотека трябва да оценява аналитично, а не да апроксимира на ръка. Компонент за електронни таблици, който претендира за съвместимост с Excel, трябва да възпроизвежда тези стойности до последната цифра, показвана от 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; // наблюдение
sh.Range['A2', 'A2'].Value := 100; // средно аритметично
sh.Range['A3', 'A3'].Value := 15; // стандартно отклонение
// Парсерът за формули на XLS използва ';' като разделител на аргументи.
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 на работната книга компилира и оценява ad-hoc формула спрямо текущия лист и връща Variant. Една подробност може да обърка потребителите при първия опит: парсерът за формули зад Calculate приема точка и запетая за разделител на аргументи, така че се пише =SUM(A1;B1), а не =SUM(A1,B1). Записаните в клетки формули запазват стандартната за Excel запетая. Същият модул за оценка изпълнява всяка статистическа функция по-долу, така че щом някоя от тях проработи в Calculate, останалите следват същия път
Двете функции, върху които е изградено всичко останало
Повечето кумулативни разпределения в този набор не се изчисляват чрез сумиране или интегриране на собствените им дефиниции. Те се изчисляват на базата на две специални функции: регуларизираната долна непълна гама функция, означена като P(a, x), и регуларизираната непълна бета функция, означена като Ix(a, b). Вътрешно това са помощните функции, на които разчита разпределителят, и веригата е къса. Кумулативната функция (CDF) на хи-квадрат е CDF на гама разпределението с параметър на формата df/2 и мащаб 2. CDF на гама разпределението е директно P(a, x). Кумулативните функции на t, F и биномното разпределение са стойности на регуларизираната непълна бета функция при съответните аргументи. CDF на Поасон е горната непълна гама функция Q. Имплементирайте добре гама и бета функциите и десетки разпределения ще наследят тяхната точност безплатно
Думата „регуларизирана“ е ключова. Суровата непълна гама функция расте като факториел, а суровият бета интеграл може да предизвика препълване (overflow) или недопълване (underflow) много преди самият отговор да го направи. Регуларизираните форми се разделят на пълната гама или бета функция, така че те живеят изцяло в интервала от нула до едно, което е точно диапазонът на вероятностите. Тази нормализация позволява на една и съща подпрограма да обслужва хи-квадрат с две степени на свобода и с двеста, без междинните членове да излизат извън границите на типа Double. Това също така обяснява защо не изчислявате CDF чрез събиране на дълга опашка от членове на плътността: всеки член носи собствена грешка от закръгляне, грешките се натрупват с изпълнението на реда, а регуларизираната специална функция избягва сумата изцяло, като оценява бързо сближаващ се ред или верижна дроб
Степенни редове под диагонала, верижни дроби над него
Подпрограмата за непълна гама функция взема едно решение преди изчислението: сравнява x с a + 1. Тази граница не е произволна. Развитието в степенен ред на P(a, x) се сближава бързо, когато x е малко спрямо a, и бавно (в крайна сметка неизползваемо), когато x е голямо. Верижната дроб има противоположното свойство. Затова енджинът използва степенния ред за x под a + 1 и верижна дроб по метода на Lentz за x при или над a + 1, като всяко разклонение изпълнява само работата, за която е подходящо
Верижната дроб се нуждае от защита. Методът на Lentz работи чрез поддържане на текущ числител и знаменател и обръщане (инвертиране) на знаменателя на всяка стъпка; ако някое от тях се доближи до нула, операцията се проваля. Решението е малък долен праг (floor): всеки път, когато междинен член падне под приблизително 1e-30 по абсолютна стойност, той се ограничава до 1e-30, което поддържа рекурсията крайна, без да нарушава сближената стойност. Същото ограничение се появява и във верижната дроб за непълната бета функция по същата причина. Това е малка константа, която върши изключително важна работа — разликата между стабилна оценка и деление на стойност, която не се отличава от нула
Горната опашка Q(a, x) е просто 1 минус P(a, x), и така се изчислява кумулативното разклонение на Поасон: вероятността за най-много k събития със средно λ е Q(k + 1, λ). Насочването му през горната непълна гама функция, вместо сумиране на k + 1 члена на Поасон, отново е избор за оценка на един сближаващ се израз вместо натрупване на множество малки стойности
Дискретни вероятности без препълване на факториела
Дискретните разпределения носят друга опасност. Вероятностната маса на биномно разпределение включва биномен коефициент, а коефициентът за избор на 26 от 52 е огромно цяло число. Ако го изчислите директно, числителят ще препълни типа Double преди делението, което би трябвало да го върне към разумна вероятност. Енджинът никога не го изчислява директно. Той изчислява факториелите в логаритмично пространство чрез логаритмичната гама функция (log-gamma), събира и изважда логаритмите, добавя логаритъма на вероятностите за успех и неуспех и накрая експоненцира веднъж най-накрая
// Вероятностна маса на биномно разпределение, изчислена изцяло в логаритмично пространство.
// LnGammaF(n+1) е ln(n!); трите логаритмични факториела съставят ln(C(n,k)),
// и целия експонент се изгражда преди едно извикване на Exp.
result := Exp(LnGammaF(nt + 1) - LnGammaF(kk + 1) - LnGammaF(nt - kk + 1)
+ kk * Ln(pp) + (nt - kk) * Ln(1 - pp));Самата логаритмична гама функция е апроксимация на Lanczos, точна по цялата положителна ос и евтина за изчисляване. Тъй като всяко голямо количество се пази като логаритъм до финалното извикване на Exp, най-голямото число, което подпрограмата създава, е самата вероятност (която е най-много едно). Вероятностната функция на Поасон следва същата рецепта, като единичният log-gamma член замества факториела в знаменателя. Граничните случаи са обособени, когато p е точно нула или едно, така че кодът никога не извиква Ln(0). HotXLS връща 0.2460938 за BINOM.DIST(5,10,0.5,FALSE) и 0.6766764 за кумулативната POISSON.DIST(2,2,TRUE), съвпадайки с Excel до последната цифра
Обратни функции чрез стесняване на интервала (bracketing) на правата крива
Само една обратна функция в този набор има бърза директна формула. NORM.S.INV, обратната стандартна нормална функция, използва рационална апроксимация на Acklam — двойка полиномни отношения, точни приблизително до прецизността на типа Double в целия диапазон, разделени на централна област и две опашки. Това е изчисление в аналитичен вид без итерация
Другите обратни функции нямат такава формула, затова енджинът ги инвертира числено. Той обгражда отговора с долна и горна граница, избрани от носителя на разпределението, и след това разделя интервала на две (бисекция): оценява правата CDF в средната точка, премества границата, която запазва целевата вероятност вътре, и повтаря до стесняване на интервала. За обратните функции на гама и хи-квадрат интервалът започва от нула и щедра горна оценка на базата на формата и мащаба, като удвоява горната граница, ако вероятността все още не е обхваната. Обратната функция на t обхваща симетрични граници, които се разширяват навън; обратната функция на F разделя неотрицателен интервал. Времевите разходи са няколко десетки изчисления на CDF на извикване, което е незабележимо за скоростта на електронна таблица, а предимството е, че всяка обратна функция е точно толкова точна, колкото и правата функция, която инвертира. Ето защо процес като CHISQ.DIST(CHISQ.INV(0.7,5),5,TRUE) връща точно 0.7
Десетичният логаритъм, който се криеше в опашката
Подпрограмата на Acklam за обратна нормална функция има три разклонения. Широкото централно разклонение, използвано при вероятности между 0.025 и 0.975, прекарва входа през полиномно отношение без никакви логаритми в него. Двете опашни разклонения за много малки или много големи вероятности първо изчисляват логаритъм от входа, тъй като опашката се държи като квадратен корен от минус естествения логаритъм от p
Ранна версия на опашното разклонение изчисляваше десетичен логаритъм (base-10 log), където трябваше да бъде естественият логаритъм. Двата се различават с константен коефициент от около 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-разпределението на Student съдържа тънкост, която лесно може да се обърка. Стойността ѝ идва от регуларизираната непълна бета функция, оценена при df / (df + x²), но тази бета стойност е вероятността в опашката отвъд мащаба на x, а не кумулативната вероятност до x. Симетричната форма на t-разпределението означава, че преобразуването зависи от това от коя страна на нулата попада x
// CDF на t-разпределението. ib е регуларизираната непълна бета при df/(df+x*x),
// която измерва симетричната опашка. Кумулативната стойност зависи от
// знака на x; връщането на ib без преобразуване дава грешна опашка.
ib := BetaIF(df / 2, 0.5, df / (df + x * x));
if x > 0 then
result := 1 - 0.5 * ib // над средното: едно минус половината от опашката
else if x < 0 then
result := 0.5 * ib // под средното: половината от опашката
else
result := 0.5; // точно при среднотоЗа x над нулата кумулативната вероятност е едно минус половината от симетричната опашка; за x под нулата е половината от тази опашка; при нула е точно една втора. Върнете ли бета стойността директно, ще отчетете грешната страна на разпределението, с разлика в цялото тяло на кривата за всяко ненулево x. Вариантите с дясна опашка и с две опашки стъпват на същото разклонение, поради което T.DIST.2T(1,1) връща 0.5, а T.DIST(1,1,TRUE) — 0.75, а обратната функция T.INV извършва бисекция спрямо тази коригирана CDF, така че процесът се затваря успешно
Нищо от това не се вижда от клетката, и това е търсеният резултат. Пишете формула и четете число, което съвпада с Excel. Ако разширявате енджина със собствена логика, механиката на регистриране на функция е обяснена в нашия преглед на формулния енджин и потребителските функции, а начинът, по который формулите достигат до други листове и именовани диапазони, е покрит в статията за именовани диапазони и препратки между листове. Всичко това се доставя в HotXLS spreadsheet component за Delphi и C++Builder, заедно с API за четене, запис, диаграми и форматиране, описани другаде в този блог