Digite =NORM.DIST(115,100,15,TRUE) em uma célula e o Excel retornará 0,8413447 sem cerimônia. A chamada parece uma busca em tabela. Mas não é. Por trás desse único número está a distribuição normal cumulativa, uma integral sem forma fechada, e por trás de CHISQ.INV.RT e BETA.DIST estão funções especiais que uma biblioteca cuidadosa precisa avaliar, não aproximar manualmente. Um componente de planilha que alega compatibilidade com o Excel deve reproduzir esses valores até o último dígito que o Excel mostra, o que significa reproduzir os métodos numéricos, e não apenas os nomes das funções.
O HotXLS implementa mais de cinquenta dessas funções estatísticas, e o trabalho que as torna corretas é quase totalmente invisível a partir da barra de fórmulas. Esta é uma visão de como o motor as calcula: o núcleo de função especial compartilhado, as decisões de ramificação que mantêm a aritmética estável e um bug de normal inversa que se escondeu na cauda por muito tempo porque o caso comum nunca tocava a linha problemática.
Uma chamada de planilha, cinquenta distribuições por trás dela
As funções abrangem as famílias que uma pasta de trabalho de estatística utiliza. Há a família normal, NORM.DIST e NORM.S.DIST com suas 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 auxiliares de inferência como CONFIDENCE.T e CONFIDENCE.NORM. Do ponto de vista de quem chama, cada uma é uma única fórmula. Você define as entradas nas células, solicita a avaliação da pasta de trabalho 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 na pasta de trabalho compila e avalia uma fórmula ad-hoc contra a planilha ativa e retorna um Variant. Um detalhe costuma atrapalhar na primeira tentativa: o analisador de fórmulas por trás de Calculate usa o ponto e vírgula como separador de argumentos, portanto é =SUM(A1;B1), não =SUM(A1,B1). As fórmulas de células armazenadas mantêm a vírgula padrão do Excel. O mesmo avaliador despacha cada função estatística abaixo, de modo que, uma vez que uma delas funcione no Calculate, as demais seguem o mesmo caminho.
As duas funções sobre as quais tudo o mais é construído
A maioria das distribuições cumulativas neste conjunto não é calculada somando ou integrando suas próprias definições. Elas são calculadas a partir de duas funções especiais: a gama incompleta inferior regularizada, escrita como P(a, x), e a beta incompleta regularizada, escrita como Ix(a, b). Internamente, estes são os auxiliares nos quais os despachantes se apoiam, e a cadeia é curta. A CDF qui-quadrado é a CDF gama com forma df/2 e escala 2. A CDF gama é P(a, x) diretamente. As funções cumulativas t, F e binomial são todas valores da beta incompleta regularizada com os argumentos corretos. A CDF de Poisson é a gama incompleta superior Q. Implemente bem as funções gama e beta e uma dúzia de distribuições herdará sua precisão gratuitamente.
A palavra "regularizada" é o ponto principal. A gama incompleta bruta cresce como um fatorial e a integral beta bruta pode sofrer subfluxo (underflow) ou estouro (overflow) muito antes do resultado final. As formas regularizadas são divididas pela gama ou beta completa, de modo que vivem inteiramente no intervalo de zero a um, que é exatamente a faixa que uma probabilidade ocupa. Essa normalização é o que permite que a mesma rotina atenda a um qui-quadrado com dois graus de liberdade e a outro com duzentos, sem que os termos intermediários ultrapassem o limite de um double. Isso também explica por que você não calcula uma CDF somando uma longa cauda de termos de densidade: cada termo carrega seu próprio erro de arredondamento, os erros se acumulam à medida que a série é executada e a função especial regularizada evita a soma por completo, avaliando em seu lugar uma série de convergência rápida ou 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 calcular qualquer coisa: ela compara x com a + 1. Esse 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, tornando-se inútil, quando x é grande. A fração contínua tem a característica oposta. Portanto, o motor usa a série de potências para x abaixo de a + 1 e uma fração contínua de Lentz para x igual ou acima de a + 1, de modo que cada ramificação execute apenas o trabalho no qual é eficiente.
A fração contínua precisa 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 e, se qualquer um deles se aproximar de zero, a inversão falha. A correção é um piso minúsculo: sempre que um termo intermediário cai abaixo de aproximadamente 1e-30 em magnitude, ele é limitado a 1e-30, o que mantém a recorrência finita sem perturbar o valor convergido. O mesmo limite aparece na fração contínua da beta incompleta pelo mesmo motivo. É uma pequena constante fazendo um trabalho estrutural importante, a diferença entre uma avaliação estável e uma divisão por algo indistinguível de zero.
A cauda superior, Q(a, x), é simplesmente 1 menos P(a, x), e é assim que a ramificação cumulativa de Poisson é calculada: a probabilidade de no máximo k eventos com média λ é Q(k + 1, λ). Direcioná-la pela gama incompleta superior em vez de somar k + 1 termos de Poisson é, novamente, uma escolha para avaliar uma única expressão convergente em vez de acumular muitas pequenas.
Massas discretas sem estouro de fatorial
As distribuições discretas apresentam um risco diferente. Uma massa de probabilidade binomial envolve um coeficiente binomial, e o coeficiente para cinquenta e dois escolhe vinte e seis é um inteiro enorme. Forme-o diretamente e o numerador estourará um double antes da divisão que o traria de volta a uma probabilidade sensata. O motor nunca o forma diretamente. Ele calcula os fatoriais em espaço logarítmico por meio da função log-gama, soma e subtrai os logs, incorpora o log das probabilidades de sucesso e fracasso e aplica a exponencial uma única vez no final.
// 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, precisa em todo o eixo positivo e barata de avaliar. Como toda quantidade grande é mantida como seu logaritmo até o Exp final, o maior número que a rotina chega a materializar é a própria probabilidade, que é no máximo um. A função de massa de Poisson segue a mesma receita, com o termo log-gama único substituindo 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 chame Ln(0). O HotXLS retorna 0,2460938 para BINOM.DIST(5,10,0.5,FALSE) e 0,6766764 para o cumulativo POISSON.DIST(2,2,TRUE), correspondendo ao Excel em todos os dígitos impressos.
Inversas por enquadramento da curva direta
Uma função de distribuição inversa faz a pergunta 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. NORM.S.INV, a normal padrão inversa, usa uma aproximação racional de Acklam, um par de razões polinomiais precisas até aproximadamente a precisão de um double em todo o intervalo, dividido em uma região central e duas caudas. É uma avaliação em forma fechada sem iteração.
As outras inversas não possuem tal fórmula, então o motor as inverte numericamente. Ele enquadra a resposta com um limite inferior e superior escolhido a partir do suporte da distribuição e, em seguida, faz a bisseção: avalia a CDF direta no ponto médio, move o limite que mantém a probabilidade alvo contida e repete até que o intervalo seja estreito. Para as inversas gama e qui-quadrado, o enquadramento começa em zero e uma estimativa superior generosa construída a partir da forma e escala, dobrando o limite superior se a probabilidade ainda não estiver contida. A inversa t enquadra limites simétricos que se alargam para fora; a inversa F faz a bisseção em um intervalo não negativo. O custo é de algumas dezenas de avaliações de CDF por chamada, o que é invisível na velocidade da planilha, e o benefício é que cada inversa é exatamente tão precisa quanto a função direta que ela inverte. É por isso que uma viagem de ida e volta como CHISQ.DIST(CHISQ.INV(0.7,5),5,TRUE) retorna 0,7 com altíssima precisão.
O logaritmo de base 10 que se escondeu na cauda
Aqui está o bug que vale a pena contar, porque é do tipo que sobrevive por muito tempo. A rotina de normal inversa de Acklam possui três ramificações. A ramificação central ampla, usada sempre que a probabilidade fica entre cerca de 0,025 e 0,975, executa a entrada por meio de uma razão polinomial sem logaritmo em lugar nenhum. As duas ramificações de cauda, para probabilidades muito pequenas ou muito grandes, tiram um logaritmo da entrada primeiro, porque a cauda se comporta como a raiz quadrada de menos o logaritmo natural de p.
Uma versão inicial da ramificação de cauda tirava um logaritmo de base 10 onde o logaritmo natural deveria estar. Os dois diferem por um fator constante de cerca de 2,30, de modo que os resultados de cauda estavam errados por uma margem consistente e considerável. E, no entanto, a função parecia boa em todas as verificações casuais, porque as verificações casuais ocorrem no meio. NORM.S.INV(0.5) é zero, NORM.S.INV(0.975) é o clássico 1,959964, e ambos passam pelo polinômio central que não chama logaritmo algum. O erro só aparecia quando uma probabilidade entrava na cauda, como NORM.S.INV(0.001), que deve retornar -3,0902323 e, em vez disso, retornava com desvio devido à proporção entre logaritmo natural e base 10. Qualquer função que dependesse da normal inversa em sua cauda, incluindo os auxiliares de intervalo de confiança, herdava o mesmo desvio. A lição é simples e custosa: uma função com estrutura de ramificação precisa de pontos de teste dentro de cada ramificação, porque um caminho comum correto ocultará alegremente um caminho raro quebrado. A correção foi a mudança de um único token do log de base 10 para o log natural, e os valores de cauda alinharam-se aos do Excel.
O sinal de x decide a cauda da distribuição t
A função cumulativa t de Student traz uma sutileza fácil de errar. Seu valor vem da beta incompleta regularizada avaliada em df / (df + x²), mas esse valor beta é a probabilidade na cauda além da magnitude de x, não a probabilidade cumulativa até x. A forma simétrica da distribuição t significa que a conversão depende de qual lado de zero x cai.
// 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 meanNada disso é visível a partir da célula, e esse é o resultado pretendido. Você escreve uma fórmula e lê um número que concorda com o Excel. Se você estiver estendendo o motor com sua própria lógica, os mecanismos de registro de uma função são abordados em nosso passo a passo sobre o motor de fórmulas e funções personalizadas, e a maneira como as fórmulas alcançam diferentes planilhas e intervalos nomeados é descrita no artigo sobre nomes definidos e fórmulas entre planilhas. Tudo isso é enviado dentro do componente de planilha HotXLS para Delphi e C++Builder, junto com as APIs de leitura, gravação, gráfico e formatação abordadas em outras partes deste blog.