Funciones estadísticas de Excel en Delphi: NORM, CHISQ, BETA

Escriba =NORM.DIST(115,100,15,TRUE) en una celda y Excel devolverá 0.8413447 sin rodeos. La llamada parece una simple búsqueda en una tabla, pero no lo es. Detrás de ese número está la distribución normal acumulada, una integral que carece de forma cerrada; y detrás de CHISQ.INV.RT y BETA.DIST residen funciones especiales que una biblioteca minuciosa debe evaluar, no aproximar a mano. Un componente de hoja de cálculo que declare compatibilidad con Excel debe reproducir estos valores hasta el último dígito que muestra Excel, lo que implica reproducir los métodos numéricos, no solo los nombres de las funciones.

HotXLS implementa más de cincuenta de estas funciones estadísticas, y el trabajo que las hace correctas es casi por completo invisible desde la barra de fórmulas. Este es un recorrido por la forma en que el motor las calcula: el núcleo compartido de funciones especiales, las decisiones de ramificación que mantienen estable la aritmética y un bug de normal inversa que permaneció oculto en la cola durante mucho tiempo porque el caso común nunca tocaba la línea con errores.

Una llamada de hoja de trabajo, cincuenta distribuciones detrás

Las funciones abarcan las familias a las que recurre un libro de trabajo de estadística. Está la familia normal, NORM.DIST y NORM.S.DIST con sus inversas; la familia gamma y chi-cuadrada, GAMMA.DIST, CHISQ.DIST, CHISQ.DIST.RT, CHISQ.INV.RT; la familia beta, BETA.DIST y BETA.INV; las distribuciones de muestreo T.DIST, T.DIST.2T, F.DIST y F.INV; el par discreto BINOM.DIST y POISSON.DIST; y los ayudantes de inferencia como CONFIDENCE.T y CONFIDENCE.NORM. Desde la perspectiva del llamador, cada una es una única fórmula: se establecen las entradas en las celdas, se solicita al libro de trabajo que las evalúe y se lee el 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;

El método Calculate del libro de trabajo compila y evalúa una fórmula ad-hoc contra la hoja activa y devuelve un Variant. Un detalle suele causar confusión en el primer intento: el analizador de fórmulas detrás de Calculate utiliza el punto y coma como separador de argumentos, por lo que es =SUM(A1;B1), no =SUM(A1,B1). Las fórmulas almacenadas en las celdas mantienen la coma estándar de Excel. El mismo evaluador despacha cada función estadística a continuación, por lo que una vez que una de ellas funciona en Calculate, las demás siguen la misma ruta.

Las dos funciones sobre las cuales se construye todo lo demás

La mayoría de las distribuciones acumuladas en este grupo no se calculan sumando o integrando sus propias definiciones. Se calculan a partir de dos funciones especiales: la función gamma incompleta inferior regularizada, denotada P(a, x), y la función beta incompleta regularizada, denotada I_x(a, b). Internamente, estos son los ayudantes en los que se apoyan los despachadores, y la cadena es corta. La CDF de chi-cuadrada es la CDF de gamma con forma df/2 y escala 2. La CDF de gamma es directamente P(a, x). Las funciones acumuladas binomial, t y F son todas valores de la función beta incompleta regularizada con los argumentos adecuados. La CDF de Poisson es la función gamma incompleta superior Q. Implemente de forma correcta las funciones gamma y beta y una docena de distribuciones heredarán su precisión de manera gratuita.

La palabra \"regularizada\" es la clave. La función gamma incompleta cruda crece como un factorial y la integral beta cruda puede sufrir desbordamiento negativo o positivo mucho antes de que lo haga el resultado. Las formas regularizadas se dividen entre la función gamma o beta completa, por lo que residen completamente en el intervalo de cero a uno, que es exactamente el rango que ocupa una probabilidad. Esa normalización es lo que permite que la misma rutina sirva para una chi-cuadrada con dos grados de libertad y para otra con doscientos sin que los términos intermedios excedan los límites de un double. También explica por qué no se calcula una CDF sumando una larga cola de términos de densidad: cada término lleva su propio error de redondeo, los errores se acumulan a medida que se ejecuta la serie, y la función especial regularizada evita la suma por completo evaluando en su lugar una serie de rápida convergencia o una fracción continua.

Serie por debajo de la diagonal, fracción continua por encima

La rutina de la función gamma incompleta toma una decisión antes de calcular cualquier cosa: compara x frente a a + 1. Ese límite no es arbitrario. La expansión en serie de potencias de P(a, x) converge rápidamente cuando x es pequeña en relación con a, y lentamente (con el tiempo de forma inútil) cuando x es grande. La fracción continua presenta el comportamiento opuesto. De modo que el motor utiliza la serie de potencias para x por debajo de a + 1 y una fracción continua de Lentz para x igual o por encima de a + 1, asignando a cada rama solo el trabajo para el que es apta.

La fracción continua necesita una protección. El método de Lentz funciona manteniendo un numerador y denominador acumulados e invirtiendo el denominador en cada paso; si alguno se acerca a cero, la inversión falla catastróficamente. La solución es un límite inferior minúsculo: cada vez que un término intermedio cae por debajo de aproximadamente 1e-30 en magnitud, se fija en 1e-30, lo que mantiene la recurrencia finita sin alterar el valor de convergencia. El mismo límite aparece en la fracción continua de la función beta incompleta por la misma razón. Es una constante pequeña que realiza un trabajo crucial: la diferencia entre una evaluación estable y una división por algo indistinguible de cero.

La cola superior, Q(a, x), es simplemente 1 menos P(a, x), y así es como se calcula la rama acumulada de Poisson: la probabilidad de obtener a lo sumo k eventos con media λ es Q(k + 1, λ). Dirigirla a través de la función gamma incompleta superior en lugar de sumar k + 1 términos de Poisson es, nuevamente, la elección de evaluar una sola expresión convergente en lugar de acumular many pequeñas.

Masas discretas sin desbordamiento factorial

Las distribuciones discretas plantean un riesgo diferente. Una masa de probabilidad binomial involucra un coeficiente binomial, y el coeficiente para 52 sobre 26 es un entero enorme. Si se calcula directamente, el numerador desborda un double antes de la división que lo devolvería a una probabilidad razonable. El motor nunca lo calcula directamente; realiza las operaciones con factoriales en el espacio logarítmico a través de la función log-gamma, suma y resta los logaritmos, incorpora el logaritmo de las probabilidades de éxito y fracaso, y calcula el valor exponencial una sola vez al 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));

La función log-gamma misma es una aproximación de Lanczos, precisa a lo largo de todo el eje positivo y de bajo costo de evaluación. Debido a que cada cantidad grande se mantiene como su logaritmo hasta la función Exp final, el número más grande que la rutina llega a materializar es la probabilidad misma, que es como máximo uno. La función de masa de Poisson sigue la misma receta, con el término único log-gamma en lugar del factorial en el denominador. Las formas cerradas se tratan de manera especial en los bordes, donde p es exactamente cero o uno, para que el código nunca llame a Ln(0). HotXLS devuelve 0.2460938 para BINOM.DIST(5,10,0.5,FALSE) y 0.6766764 para la acumulada POISSON.DIST(2,2,TRUE), coincidiendo con Excel en los dígitos mostrados.

Inversas mediante acotación de la curva directa

Una función de distribución inversa plantea la pregunta opuesta: encontrar el valor cuya CDF sea igual a ella. Solo una inversa en este grupo posee una fórmula directa rápida. NORM.S.INV, la inversa normal estándar, utiliza una aproximación racional de Acklam, un par de relaciones polinómicas precisas aproximadamente hasta la resolución de un double en todo el rango, divididas en una región central y dos colas. Es una evaluación de forma cerrada sin iteración.

Las otras inversas no tienen tal fórmula, por lo que el motor las invierte numéricamente. Acota la respuesta con un límite inferior y superior elegidos a partir del soporte de la distribución, y luego aplica bisección: evalúa la CDF directa en el punto medio, desplaza el límite que mantenga encerrada la probabilidad objetivo y repite hasta que el intervalo sea estrecho. Para las inversas de gamma y chi-cuadrada, la acotación comienza en cero y en una estimación superior generosa construida a partir de la forma y la escala, duplicando el límite superior si la probabilidad aún no está encerrada. La inversa de t acota límites simétricos que se ensanchan hacia afuera; la de F aplica bisección en un intervalo no negativo. El costo es de unas pocas docenas de evaluaciones de CDF por llamada, lo cual es imperceptible a la velocidad de una hoja de cálculo, y el beneficio es que cada inversa es exactamente tan precisa como la función directa que invierte. Es por esto que un viaje de ida y vuelta como CHISQ.DIST(CHISQ.INV(0.7,5),5,TRUE) devuelve 0.7 con gran precisión.

El logaritmo en base 10 que se escondía en la cola

Una versión temprana de la rama de la cola tomaba un logaritmo en base 10 donde correspondía el logaritmo natural. Ambos difieren en un factor constante de aproximadamente 2.30, por lo que los resultados de la cola eran incorrectos por un margen constante y considerable. Y sin embargo, la función parecía correcta en cada verificación informal, porque las verificaciones informales se concentran en el medio. NORM.S.INV(0.5) es cero, NORM.S.INV(0.975) es el valor de manual 1.959964, y ambos pasan por el polinomio central que nunca invoca un logaritmo. El error solo aparecía una vez que una probabilidad cruzaba a una cola, por ejemplo NORM.S.INV(0.001), que debe devolver -3.0902323 y en su lugar devolvía un valor alterado por la relación entre logaritmo natural y base 10. Cualquier función que dependiera de la normal inversa en su cola, incluidos los asistentes de intervalo de confianza, heredaba el mismo sesgo. La lección es común y costosa: una función con una estructura de ramas necesita puntos de prueba dentro de cada rama, porque una ruta común correcta enmascarará alegremente una ruta inusual con fallas. La corrección consistió en cambiar un solo token del logaritmo en base 10 al logaritmo natural, y los valores de la cola se ajustaron a los de Excel.

El signo de x decide la cola de la distribución t

La función acumulada t de Student presenta una sutileza fácil de calcular al revés. Su valor proviene de la función beta incompleta regularizada evaluada en df / (df + x²), pero ese valor beta es la probabilidad en la cola más allá de la magnitud de x, no la probabilidad acumulada hasta x. La forma simétrica de la distribución t implica que la conversión depende de en qué lado del cero se ubique 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

Para x por encima de cero, la probabilidad acumulada es uno menos la mitad de la cola simétrica; para x por debajo de cero, es la mitad de esa cola; en cero, es exactamente un medio. Si devuelve el valor beta directamente, informará el lado incorrecto de la distribución, fallando por todo el cuerpo de la curva para cualquier x distinto de cero. Las variantes de cola derecha y de dos colas se construyen sobre la misma rama, razón por la cual T.DIST.2T(1,1) devuelve 0.5 y T.DIST(1,1,TRUE) devuelve 0.75, y la inversa T.INV aplica bisección contra esta CDF corregida para completar el ciclo.

Nada de esto es visible desde la celda, y ese es el resultado previsto. Usted escribe una fórmula y lee un número que coincide con Excel. Si está extendiendo el motor con su propia lógica, los mecanismos para registrar una función se tratan en nuestro recorrido por el motor de fórmulas y funciones personalizadas, y la forma en que las fórmulas acceden a otras hojas y rangos nombrados se cubre en el artículo sobre nombres definidos y fórmulas entre hojas. Todo esto se incluye dentro del componente de hoja de cálculo HotXLS para Delphi y C++Builder, junto con las API de lectura, escritura, gráficos y formato tratadas en otras secciones de este blog.