Digite =NORM.DIST(115,100,15,TRUE) numa célula e o Excel devolve 0.8413447 sem demoras. A chamada parece uma simples pesquisa (lookup), mas não é. Por trás desse número reside a distribuição normal cumulativa, um integral sem forma fechada, e por trás de CHISQ.INV.RT e BETA.DIST encontram-se funções especiais que uma biblioteca rigorosa tem de avaliar, e não aproximar manualmente. Um componente de folha de cálculo que reivindique compatibilidade com o Excel tem de reproduzir estes valores até ao último dígito apresentado pelo Excel, o que implica reproduzir os métodos numéricos e não apenas os nomes das funções
O HotXLS implementa mais de cinquenta destas funções estatísticas, e o trabalho que garante a sua exatidão é quase totalmente invisível a partir da barra de fórmulas. Esta é uma explicação detalhada de como o motor as calcula: o núcleo partilhado de funções especiais, as decisões de ramificação que mantêm a estabilidade aritmética e um bug na distribuição normal inversa que permaneceu oculto na cauda (tail) durante muito tempo porque os cenários comuns nunca acediam à linha com erro
Uma chamada na folha de cálculo, vinte distribuições por trás dela
As funções abrangem as famílias habitualmente procuradas num livro de estatística. Incluem a família normal, com NORM.DIST e NORM.S.DIST e respetivas inversas; a família gama e qui-quadrado, com GAMMA.DIST, CHISQ.DIST, CHISQ.DIST.RT e CHISQ.INV.RT; a família beta, com BETA.DIST e BETA.INV; as distribuições de amostragem T.DIST, T.DIST.2T, F.DIST e F.INV; o par discreto BINOM.DIST e POISSON.DIST; e os assistentes de inferência como CONFIDENCE.T e CONFIDENCE.NORM. Do ponto de vista do programador, cada uma constitui uma única fórmula: define as entradas nas células, solicita a avaliação ao livro e lê o resultado
O método Calculate no livro compila e avalia uma fórmula ad-hoc contra a folha ativa e devolve um tipo Variant. Um detalhe costuma surpreender os usuários na primeira tentativa: o analisador de fórmulas por trás de Calculate adota o ponto e vírgula como separador de argumentos, pelo que se escreve =SUM(A1;B1) e não =SUM(A1,B1). As fórmulas guardadas nas células mantêm a vírgula padrão do Excel. O mesmo avaliador processa todas as funções estatísticas listadas abaixo; uma vez que uma destas funcione em Calculate, as restantes seguem o mesmo caminho
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;As duas funções sobre as quais tudo o resto é construído
A maioria das distribuições cumulativas neste conjunto não são calculadas somando ou integrando as suas próprias definições. São calculadas a partir de duas funções especiais: a função gama incompleta inferior regularizada, escrita como P(a, x), e a função beta incompleta regularizada, escrita como Ix(a, b). Internamente, estes são os assistentes nos quais os despachantes se apoiam, e a cadeia é curta: a CDF (função de distribuição cumulativa) qui-quadrado é a CDF gama com forma df/2 e escala 2; a CDF gama é diretamente P(a, x); as funções cumulativas t, F e binomial são todas valores da função beta incompleta regularizada com os argumentos adequados; e a CDF de Poisson é a função gama incompleta superior Q. Implemente as funções gama e beta corretamente e uma dúzia de distribuições herdará a sua exatidão gratuitamente
A palavra "regularizada" é o aspeto chave. A função gama incompleta pura cresce como um fatorial e o integral beta puro pode sofrer underflow ou overflow muito antes de o resultado ser alcançado. As formas regularizadas são divididas pela função gama ou beta completa, pelo que se situam inteiramente no intervalo de zero a um, que corresponde exatamente ao intervalo que uma probabilidade ocupa. Essa normalização é o que permite à mesma rotina processar um qui-quadrado com dois graus de liberdade e um com duzentos sem que os termos intermédios ultrapassem o limite de um double. Explica também por que não se calcula uma CDF somando uma longa cauda de termos de densidade: cada termo carrega o seu próprio erro de arredondamento, os erros acumulam-se à medida que a série avança, e a função especial regularizada contorna a soma por completo, avaliando em vez disso uma série de convergência rápida ou uma fração contínua
Série abaixo da diagonal, fração contínua acima dela
A rotina gama incompleta toma uma decisão antes de efetuar qualquer cálculo: compara x com a + 1. Este limite não é arbitrário: a expansão em série de potências de P(a, x) converge rapidamente quando x é pequeno em relação a a, e lentamente (acabando por ser inútil) quando x é grande. A fração contínua possui o comportamento oposto. Desta forma, o motor utiliza a série de potências para valores de x abaixo de a + 1 e a fração contínua de Lentz para x igual ou superior a a + 1, garantindo que cada ramificação execute apenas o trabalho no qual é eficiente
A fração contínua necessita de uma proteção: o método de Lentz funciona mantendo um numerador e denominador dinâmicos e invertendo o denominador a cada passo; se algum deles se aproximar de zero, a inversão falha. A solução consiste num limite mínimo minúsculo: sempre que o módulo de um termo intermédio cai abaixo de cerca de 1e-30, é fixado em 1e-30, mantendo a recorrência finita sem perturbar o valor de convergência. O mesmo limite surge na fração contínua da função beta incompleta pelo mesmo motivo. Trata-se de uma pequena constante a realizar um trabalho fundamental, ditando a diferença entre uma avaliação estável e uma divisão por um valor indistinguível de zero
A cauda superior, Q(a, x), corresponde simplesmente a 1 menos P(a, x), e é desta forma que a ramificação cumulativa de Poisson é calculada: a probabilidade de ocorrer no máximo k eventos com média λ é Q(k + 1, λ). Encaminhar esta operação através da gama incompleta superior em vez de somar k + 1 termos de Poisson constitui, mais uma vez, a opção de avaliar uma única expressão convergente em vez de acumular múltiplos valores pequenos
Massas discretas sem overflow de fatorial
As distribuições discretas trazem um perigo diferente. Uma massa de probabilidade binomial envolve um coeficiente binomial, e o coeficiente para cinquenta e dois combinados vinte a vinte e seis é um número inteiro gigante. Se o gerar diretamente, o numerador sofrerá overflow num double antes da divisão que o traria de volta a uma probabilidade lógica. O motor nunca o gera dessa forma: calcula os fatoriais em espaço logarítmico através da função log-gama, adiciona e subtrai os logaritmos, junta o logaritmo das probabilidades de sucesso e falha e aplica a exponencial uma única vez no final
A própria função log-gama é uma aproximação de Lanczos, exata ao longo de todo o eixo positivo e computacionalmente económica. Como cada quantidade grande é mantida como o seu logaritmo até à chamada final de Exp, o maior número que a rotina chega a materializar é a própria probabilidade, cujo valor máximo é um. A função de massa de Poisson segue o mesmo método, com o termo log-gama a substituir o fatorial no denominador. Os casos limites são tratados de forma especial nas extremidades, onde p é exatamente zero ou um, garantindo que o código nunca chame Ln(0). O HotXLS devolve 0.2460938 para BINOM.DIST(5,10,0.5,FALSE) e 0.6766764 para a cumulativa POISSON.DIST(2,2,TRUE), correspondendo aos dígitos apresentados pelo Excel
// 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));Inversas por enquadramento da curva direta
Uma função de distribuição inversa coloca a questão oposta: dada uma probabilidade, encontre o valor cuja CDF seja igual a ela. Apenas uma inversa neste conjunto possui uma fórmula direta rápida: a NORM.S.INV (normal padrão inversa) recorre a uma aproximação racional de Acklam, um par de razões polinomiais exatas até perto da precisão de um double em todo o intervalo, dividida numa região central e duas caudas. Trata-se de uma avaliação de forma fechada, sem iterações
As restantes inversas não dispõem de tal fórmula, pelo que o motor as inverte numericamente. Enquadra o resultado com um limite inferior e superior escolhidos a partir do suporte da distribuição e, em seguida, aplica a bisseção: avalia a CDF direta no ponto médio, move o limite que mantém a probabilidade alvo enquadrada e repete o processo até o intervalo ser estreito. Para as inversas gama e qui-quadrado, o enquadramento começa em zero e numa estimativa superior generosa construída com base na forma e escala, duplicando o limite superior se a probabilidade ainda não estiver enquadrada. A inversa t enquadra limites simétricos que se alargam para o exterior; a inversa F realiza a bisseção num intervalo não negativo. O custo reside em algumas dezenas de avaliações de CDF por chamada (o que é impercetível à velocidade da folha de cálculo), e a vantagem é que cada inversa é exatamente tão exata quanto a função direta que inverte. É por isso que um percurso de ida e volta como CHISQ.DIST(CHISQ.INV(0.7,5),5,TRUE) devolve 0.7 com elevada precisão
O logaritmo de base 10 que se escondeu na cauda
A rotina de normal inversa de Acklam possui três ramificações: a ramificação central ampla, utilizada sempre que a probabilidade se situa entre cerca de 0.025 e 0.975, processa a entrada através de uma razão polinomial sem qualquer logaritmo; as duas ramificações de cauda, para probabilidades muito pequenas ou muito grandes, extraem primeiro um logaritmo da entrada, visto que o comportamento da cauda se assemelha à raiz quadrada do simétrico do logaritmo natural de p
Uma versão inicial da ramificação de cauda adotava um logaritmo de base 10 onde devia estar o logaritmo natural. Os dois diferem por um fator constante de cerca de 2.30, pelo que os resultados na cauda apresentavam um erro consistente e significativo. Contudo, a função parecia correta em qualquer verificação superficial, uma vez que estas se focavam na zona central: NORM.S.INV(0.5) é zero, NORM.S.INV(0.975) corresponde ao valor clássico 1.959964, e ambos os casos correm no polinómio central que nunca invoca logaritmos. O erro apenas surgia quando a probabilidade entrava na cauda, por exemplo em NORM.S.INV(0.001), que deve devolver -3.0902323 e, em vez disso, apresentava um desvio proporcional à razão entre o logaritmo natural e o de base 10. Qualquer função dependente da normal inversa na sua cauda, incluindo os assistentes de intervalo de confiança, herdavam a mesma distorção. A lição é simples mas dispendiosa: uma função com estrutura de ramificação necessita de pontos de teste em cada ramo, visto que um caminho comum correto ocultará com facilidade um ramo raramente acedido com erros. A correção consistiu na alteração de um único token do logaritmo de base 10 para o logaritmo natural, fazendo os valores na cauda corresponderem aos do Excel
O sinal de x determina a cauda da distribuição t
A função cumulativa t de Student apresenta uma subtileza fácil de inverter: o seu valor provém da função beta incompleta regularizada avaliada em df / (df + x²), mas esse valor beta corresponde à probabilidade na cauda além do módulo de x, e não à probabilidade cumulativa até x. A forma simétrica da distribuição t implica que a conversão dependa do lado de zero onde x se situa
Para x superior a zero, a probabilidade cumulativa é um menos metade da cauda simétrica; para x inferior a zero, corresponde a metade dessa cauda; em zero, é exatamente metade. Se devolver o valor beta diretamente, reportará o lado incorreto da distribuição, com um desvio equivalente a todo o corpo da curva para qualquer x diferente de zero. As variantes de cauda direita e de duas caudas assentam na mesma ramificação, razão pela qual T.DIST.2T(1,1) devolve 0.5 e T.DIST(1,1,TRUE) devolve 0.75, e a inversa T.INV realiza a bisseção contra esta CDF corrigida para que a correspondência de ida e volta se verifique
Nada disto é visível a partir da célula, e esse é o resultado pretendido: escreve uma fórmula e lê um número que corresponde ao do Excel. Se estiver a estender o motor com a sua própria lógica, o mecanismo de registro de funções é abordado no artigo explicativo sobre o motor de fórmulas e funções personalizadas, e a forma como as fórmulas acedem a folhas e intervalos nomeados é detalhada no artigo sobre nomes definidos e fórmulas entre folhas. Todo este conjunto de funcionalidades é fornecido no HotXLS spreadsheet component para Delphi e C++Builder, juntamente com as APIs de leitura, escrita, gráficos e formatação abordadas noutras secções deste blog
// 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