Статистически функции на Excel в Delphi: NORM, CHISQ, BETA

Въведете =NORM.DIST(115,100,15,TRUE) в клетка и Excel връща 0.8413447 без никакви формалности. Извикването изглежда как търсене в таблица. Не е. Зад това единствено число стои кумулативното нормално разпределение - интеграл без затворена форма, а зад 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;   // 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 на работната книга компилира и оценява ad-hoc формула спрямо живия лист и връща Variant. Една подробност спъва хората при първия опит: парсерът на формули зад Calculate приема точка и запетая за свой разделител на аргументи, така че се пише =SUM(A1;B1), а не =SUM(A1,B1). Съхранените формули на клетки запазват стандартната за Excel запетая. Същият оценител изпраща всяка статистическа функция по-долу, така че след как една от тях заработи в Calculate, останалите следват същия път.

Двете функции, върху които е изградено всичко останало

Повечето кумулативни разпределения в този набор не се изчисляват чрез сумиране или интегриране на собствените им дефиниции. Те се изчисляват от две специални функции: регулираната долна непълна гама функция (regularized lower incomplete gamma), изписана как P(a, x), и регулираната непълна бета функция (regularized incomplete beta), изписана как I_x(a, b). Вътрешно това са помощниците, на които се опират диспечерите, а веригата е кратка. CDF (функцията на разпределение) на хи-квадрат е гама CDF с форма df/2 и мащаб 2. Гама CDF е директно P(a, x). Кумулативните функции на t, F и биномното разпределение са стойности на регулираната непълна бета при правилните аргументи. CDF на Поасон е горната непълна гама Q. Имплементирайте добре гама и бета функциите и дузина разпределения ще наследят безплатно тяхната точност.

Думата "регулирана" (regularized) е цялата същност. Суровата непълна гама расте как факториел, а суровият бета интеграл може да се провали от недопълване (underflow) или препълване (overflow) дълго преди отговорът да го направи. Регулираните форми са разделени на пълната гама или бета, така че те живеят изцяло в интервала от нула до едно, което е точно диапазонът, който заема една вероятност. Тази нормализация е това, което позволява на една и съща рутина да обслужва хи-квадрат с две степени на свобода и такава с двеста, без междинните членове да излизат извън границите на double. Тя също така обяснява защо не изчислявате CDF чрез събиране на дълга опашка от членове на плътността: всеки член носи собствена грешка от закръгляне, грешките се натрупват с изпълнението на реда, а регулираната специална функция избягва изцяло сумата, като вместо това оценява бързо сходящ ред или верижна дроб.

Редове под диагонала, верижни дроби над него

Подпрограмата за непълна гама взема едно решение, преди да изчисли каквото и да било: тя сравнява x с a + 1. Тази граница не е произволна. Развитието в степенен ред на P(a, x) се сближава бързо, когато x е малко спрямо a, и бавно (в крайна сметка безполезно), когато x е голямо. Верижната дроб има обратния характер. Така че машината използва степенния ред за x под a + 1 и верижната дроб на Lentz за x на или над a + 1, като всеки клон е помолен да върши само работата, в която е добър.

Верижната дроб се нуждае от една защита. Методът на Lentz работи чрез носене на текущ числител и знаменател и обръщане на знаменателя на всяка стъпка, и ако някое от тях се доближи до нула, обръщането се проваля. Решението е малък долен праг: винаги когато междинен член падне под приблизително 1e-30 по величина, той се фиксира на 1e-30, което поддържа рекурентната връзка крайна, без да нарушава сближената стойност. Същото фиксиране се появява и във верижната дроб на непълната бета по същата причина. Това е малка константа, вършеща важна носеща работа - разликата между стабилна оценка и деление на нещо, неразличимо от нула.

Горната опашка, Q(a, x), е просто 1 минус P(a, x), и по този начин се изчислява кумулативният клон на Поасон: вероятността за най-много k събития със средно λ е Q(k + 1, λ). Насочването му през горната непълна гама, вместо сумирането на k + 1 членове на Поасон, отново е избор да се оцени един сходящ израз, вместо да се натрупват много малки такива.

Дискретни маси без препълване на факториели

Дискретните разпределения пораждат различна опасност. Биномната вероятностна маса включва биномен коефициент, а коефициентът за "петдесет и две над двадесет и шест" е огромно цяло число. Създайте го директно и числителят препълва тип double преди делението, което би го върнало обратно до разумна вероятност. Машината никога не го изгражда директно. Тя изчислява факториелите в логаритмично пространство чрез логаритмичната гама функция (log-gamma function), събира и изважда логаритмите, добавя логаритъма на вероятностите за успех и неуспех и извършва експоненциране точно веднъж накрая.

// 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));

Самата логаритмична гама функция е апроксимация на Ланцош (Lanczos approximation), точна по цялата положителна ос и евтина за оценка. Тъй като всяко голямо количество се държи как свой логаритъм до крайното извикване на Exp, най-голямото число, което рутината някога материализира, е самата вероятност, която е най-много едно. Масовата функция на Поасон следва същата рецепта, като единичният член за логаритмична гама замества факториела в знаменателя. Затворените форми се разглеждат как специални случаи по краищата, където p е точно нула или едно, така че кодът никога не извиква Ln(0). HotXLS връща 0.2460938 за BINOM.DIST(5,10,0.5,FALSE) и 0.6766764 за кумулативното POISSON.DIST(2,2,TRUE), съвпадайки с Excel във всички изписани цифри.

Инверсии чрез заграждане (bracketing) на правата крива

Инверсната функция на разпределение задава противоположния въпрос: при дадена вероятност, намерете стойността, чиято кумулативна функция (CDF) е равна на нея. Само една инверсия в този набор има бърза директна формула. NORM.S.INV, инверсното стандартно нормално разпределение, използва рационална апроксимация на Acklam - двойка полиномни съотношения, точни до приблизително прецизността на тип double в целия диапазон, разделени на централна област и две опашки. Това е оценка в затворена форма без итерация.

Другите инверсии нямат такава формула, така че машината ги обръща числено. Тя загражда (brackets) отговора с долна и горна граница, избрани от поддръжката на разпределението, и след това разделя наполовина (bisects): оценява правата кумулативна функция в средната точка, премества границата, която поддържа целевата вероятност затворена, и повтаря, докато интервалът се стесни. За инверсиите на гама и хи-квадрат заграждането започва от нула и щедра горна оценка, изградена от формата и мащаба, като удвоява горната граница, ако вероятността все още не е затворена. Инверсията на t загражда симетрични граници, които се разширяват навън; инверсията на F разделя наполовина върху неотрицателен интервал. Цената е няколко десетки оценки на кумулативната функция на извикване, което е невидимо при скоростта на електронна таблица, а ползата е, че всяка инверсия е точно толкова точна, колкото и правата функция, която обръща. Ето защо двупосочно пътуване как 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-разпределението на Студент носи тънкост, която е лесно да се обърка. Стойността му идва от регулираната непълна бета, оценена при 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 spreadsheet component за Delphi and C++Builder, заедно с API за четене, писане, диаграми и форматиране, разгледани на други места в този блог.