Digite =NORM.DIST(115,100,15,TRUE) numa célula e o Excel devolve 0.8413447 sem demoras. A chamada parece uma simples consulta. Mas não é. Por trás desse número único está 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 deve reproduzir estes valores até ao último dígito que o Excel apresenta, o que implica reproduzir os métodos numéricos, e não apenas os nomes das funções.
HotXLS implementa mais de cinquenta destas funções estatísticas, e o trabalho que garante a sua exatidão é quase por completo invisível na barra de fórmulas. Esta é uma análise de como o motor as calcula: o núcleo partilhado de funções especiais, as decisões de ramificação que mantêm a aritmética estável, e um erro na distribuição normal inversa que permaneceu oculto na cauda da curva durante muito tempo porque o caso de uso comum nunca acedia à linha com problemas.
Uma chamada de folha de cálculo, cinquenta distribuições por trás
As funções abrangem as famílias que um livro de trabalho de estatística utiliza. Existe a família normal, NORM.DIST e NORM.S.DIST com as respetivas inversas; a família gama e qui-quadrado, GAMMA.DIST, CHISQ.DIST, CHISQ.DIST.RT, CHISQ.INV.RT; a família beta, 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 de quem efetua a chamada, cada uma constitui uma única fórmula. Define as entradas nas células, solicita ao livro que as avalie e lê o resultado.
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;O método Calculate no livro compila e avalia uma fórmula ad-hoc face à folha ativa e devolve um Variant. Um pormenor costuma induzir em erro na primeira tentativa: o analisador de fórmulas por trás de Calculate aceita o ponto e vírgula como o seu separador de argumentos, pelo que se escreve =SUM(A1;B1), e não =SUM(A1,B1). As fórmulas de células guardadas mantêm a vírgula padrão do Excel. O mesmo avaliador processa todas as funções estatísticas descritas abaixo, pelo que se uma destas funcionar no Calculate, as restantes seguirão o mesmo caminho.
As duas funções nas quais tudo o resto se baseia
A maioria das distribuições cumulativas neste conjunto não é calculada somando ou integrando as suas próprias definições. São calculadas a partir de duas funções especiais: a gama incompleta inferior regularizada, representada por P(a, x), e a beta incompleta regularizada, escrita como Ix(a, b). Internamente, estes são os assistentes em que os processadores se apoiam, e a cadeia de dependências é curta. A CDF 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. A CDF de Poisson é a gama incompleta superior Q. Implemente as funções gama e beta corretamente e uma dúzia de distribuições herdará a sua exatidão sem esforço adicional.
O termo "regularizada" traduz-se no aspeto fundamental. A gama incompleta bruta cresce como um fatorial e o integral beta bruto pode sofrer transbordo negativo (underflow) ou positivo (overflow) muito antes de o resultado final se consolidar. As formas regularizadas são divididas pela gama ou beta completa, situando-se inteiramente no intervalo de zero a um, que é exatamente o espectro abrangido por uma probabilidade. Essa normalização é o que permite que a mesma rotina sirva um qui-quadrado com dois graus de liberdade e outro 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 sequência de termos de densidade: cada termo transporta o seu próprio erro de arredondamento, os erros acumulam-se à medida que a série corre, 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 de gama incompleta toma uma decisão antes de efetuar qualquer cálculo: compara x com a + 1. Esse limite não é arbitrário. A expansão por série de potências de P(a, x) converge rapidamente quando x é pequeno em relação a a, e lentamente, tornando-se eventualmente inútil, quando x é grande. A fração contínua possui a característica oposta. Assim, o motor utiliza a série de potências para x abaixo de a + 1 e uma fração contínua de Lentz para x igual ou superior a a + 1, aplicando em cada ramificação apenas o processamento para o qual ela é vocacionada.
A fração contínua necessita de uma proteção. O método de Lentz funciona mantendo um numerador e um 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 infinitesimal: sempre que um termo intermédio desce abaixo de cerca de 1e-30 em magnitude, é fixado em 1e-30, o que mantém a recorrência finita sem perturbar o valor convergente. A mesma limitação surge na fração contínua da beta incompleta pela mesma razão. Trata-se de uma constante reduzida que executa um trabalho estrutural crítico, representando a diferença entre uma avaliação estável e uma divisão por um valor indistinguível de zero.
A cauda da curva, 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 no máximo k eventos com média λ é Q(k + 1, λ). Encaminhar o cálculo através da gama incompleta superior em vez de somar k + 1 termos de Poisson traduz-se, de novo, na opção de avaliar uma única expressão convergente em vez de acumular múltiplas parcelas reduzidas.
Massas discretas sem transbordo de fatorial
As distribuições discretas apresentam um risco distinto. Uma massa de probabilidade binomial envolve um coeficiente binomial, e o coeficiente para cinquenta e dois combinado vinte e seis a vinte e seis é um inteiro gigantesco. Se for gerado diretamente, o numerador transborda a capacidade de um double antes da divisão que o traria de volta a uma probabilidade razoável. O motor nunca o gera dessa forma. Calcula os fatoriais no espaço logarítmico através da função log-gama, soma e subtrai os logaritmos, integra o logaritmo das probabilidades de sucesso e falha, e aplica a exponencial uma única vez no final de todo o processo.
// 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));A própria função log-gama é uma aproximação de Lanczos, exata ao longo de todo o eixo positivo e computacionalmente barata de avaliar. Como cada quantidade elevada é mantida como o seu logaritmo até ao Exp final, o maior número que a rotina alguma vez materializa é a própria probabilidade, que é no máximo um. A função de massa de Poisson segue a mesma lógica, com o termo único log-gama a substituir o fatorial no denominador. As formas fechadas são tratadas como casos especiais nos limites, onde p é exatamente zero ou um, para que o código nunca execute 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), matching Excel em todos os dígitos impressos.
Inversas por enquadramento da curva direta
Uma função de distribuição inversa coloca a questão oposta: dada uma probabilidade, localizar o valor cuja CDF corresponde a ela. Apenas uma inversa neste conjunto dispõe de uma fórmula direta rápida. NORM.S.INV, a distribuição normal padrão inversa, utiliza uma aproximação racional de Acklam, um par de rácios polinomiais exatos até perto da precisão de um double ao longo de todo o intervalo, dividido numa região central e duas caudas. Trata-se de uma avaliação de forma fechada sem qualquer iteração.
As outras inversas não possuem essa fórmula, pelo que o motor as inverte numericamente. Enquadra a resposta com um limite inferior e um superior escolhidos a partir do suporte da distribuição, efetuando depois uma bisseção: avalia a CDF direta no ponto médio, move o limite que mantém a probabilidade de destino enquadrada e repete o processo até o intervalo ser reduzido. Para as inversas gama e qui-quadrado, o enquadramento começa em zero e numa estimativa superior generosa baseada na forma e na escala, duplicando o limite superior se a probabilidade ainda não estiver enquadrada. A inversa t enquadra limites simétricos que se alargam para fora; a inversa F realiza a bisseção num intervalo não negativo. O custo traduz-se em algumas dezenas de avaliações de CDF por chamada, o que é impercetível na velocidade de processamento da folha de cálculo, e o benefício é que cada inversa é exatamente tão exata quanto a função direta que inverte. É por esta razão que um percurso de ida e volta como CHISQ.DIST(CHISQ.INV(0.7,5),5,TRUE) devolve 0.7 com total exatidão.
O logaritmo de base 10 que se ocultava na cauda da curva
Uma versão inicial da ramificação de cauda aplicava um logaritmo de base 10 onde devia constar o logaritmo natural. Os dois diferem por um fator constante de cerca de 2.30, pelo que os resultados na cauda da curva estavam incorretos por uma margem consistente e assinalável. No entanto, a função parecia correta em todas as verificações superficiais, porque estas incidem sobre a parte central. NORM.S.INV(0.5) é zero, NORM.S.INV(0.975) é o clássico 1.959964, e ambos correm no polinómio central que nunca recorre a logaritmos. O erro apenas surgia quando uma probabilidade entrava na cauda, como NORM.S.INV(0.001), que deve devolver -3.0902323 e, em vez disso, regressava com a diferença decorrente da relação entre o logaritmo natural e o de base 10. Qualquer função dependente da distribuição normal inversa na sua cauda, incluindo os assistentes de intervalo de confiança, herdava o mesmo desvio. A lição é simples e valiosa: uma função com estrutura de ramificações necessita de pontos de teste em cada ramificação, visto que um caminho comum correto ocultará sem dificuldade uma ramificação excecional incorreta. A correção consistiu na alteração de um único termo, passando do logaritmo de base 10 para o logaritmo natural, fazendo com que os valores da cauda passassem a coincidir com os do Excel.
O sinal de x define a cauda da distribuição t
// 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 meanPara x acima de zero, a probabilidade cumulativa corresponde a um menos metade da cauda simétrica; para x abaixo de zero, equivale a metade dessa cauda; em zero, é exatamente metade. Devolver o valor beta diretamente reportaria o lado incorreto da distribuição, desviando-se por toda a amplitude 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 efetua a bisseção face a esta CDF corrigida para que a viagem de ida e volta se conclua.
Nada disto é visível a partir da célula, e esse é o resultado pretendido. O programador escreve uma fórmula e lê um número que coincide com o Excel. Se estiver a estender o motor com a sua própria lógica, o processo de registo de uma função é coberto no nosso guia do motor de fórmulas e funções personalizadas, e o modo como as folhas acedem a outras folhas e intervalos nomeados é tratado no artigo sobre nomes definidos e fórmulas entre folhas. Tudo isto é fornecido no HotXLS spreadsheet component para Delphi e C++Builder, juntamente com as APIs de leitura, escrita, gráficos e formatação tratadas noutros locais deste blogue.