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 ceremonias. La llamada se lee como una búsqueda. No lo es. Detrás de ese número está la distribución normal acumulada, una integral sin forma cerrada, y detrás de CHISQ.INV.RT y BETA.DIST se encuentran funciones especiales que una librería cuidadosa tiene que evaluar, no aproximar a mano. Un componente de hoja de cálculo que afirme compatibilidad con Excel tiene que reproducir estos valores hasta el último dígito que muestra Excel, lo que significa 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 completamente invisible desde la barra de fórmulas. Este es un recorrido por cómo las calcula el motor: el núcleo de función especial compartido, las decisiones de rama que mantienen estable la aritmética y un error de normal inversa que se ocultó en el extremo durante mucho tiempo porque el caso común nunca tocó la línea rota.

Una llamada de hoja de cálculo, 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-cuadrado, 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 las herramientas de inferencia como CONFIDENCE.T y CONFIDENCE.NORM. Desde la posición del llamador, cada una es una sola fórmula. Establece las entradas en las celdas, solicita al libro de trabajo que evalúe y 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 en el libro de trabajo compila y evalúa una fórmula ad-hoc contra la hoja activa y devuelve un Variant. Un detalle hace tropezar a la gente en el primer intento: el analizador de fórmulas detrás de Calculate toma el punto y coma como su separador de argumentos, por lo que es =SUM(A1;B1), no =SUM(A1,B1). Las fórmulas de celdas las almacenan conservando la coma estándar de Excel. El mismo evaluador distribuye cada función estadística a continuación, por lo que una vez que una de estas funciona en Calculate, las demás siguen el mismo camino.

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

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

La palabra "regularizada" es todo el propósito. La función gamma incompleta sin procesar crece como un factorial y la integral beta sin procesar puede sufrir desbordamiento negativo o positivo mucho antes de que lo haga la respuesta. Las formas regularizadas se dividen por la función gamma o beta completa, por lo que viven 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-cuadrado con dos grados de libertad y una con doscientos sin que los términos intermedios se desborden 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 al evaluar en su lugar una serie rápidamente convergente o una fracción continua.

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

La rutina de gamma incompleta toma una decisión antes de calcular nada: compara x con 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, eventualmente de forma inútil, cuando x es grande. La fracción continua tiene el carácter opuesto. Así que el motor utiliza la serie de potencias para x por debajo de a + 1 y una fracción continua de Lentz para x en o por encima de a + 1, y a cada rama se le pide que haga solo el trabajo en el que es buena.

La fracción continua necesita una protección. El método de Lentz funciona llevando un numerador y denominador en ejecución e invirtiendo el denominador en cada paso, y si alguno se acerca a cero, la inversión explota. 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 a 1e-30, lo que mantiene la recurrencia finita sin alterar el valor convergido. La misma fijación aparece en la fracción continua de la beta incompleta por la misma razón. Es una constante pequeña que realiza un trabajo de soporte de carga, 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 acumulativa de Poisson: la probabilidad de a lo sumo k eventos con media λ es Q(k + 1, λ). Dirigirla a través de la gamma incompleta superior en lugar de sumar k + 1 términos de Poisson es, nuevamente, una opción para evaluar una expresión convergente en lugar de acumular many pequeñas.

Masas discretas sin desbordamiento factorial

Las distribuciones discretas plantean un peligro diferente. Una masa de probabilidad binomial implica un coeficiente binomial, y el coeficiente para cincuenta y dos sobre veintiséis 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 de esa forma. Calcula los 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 realiza la 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 en sí es una aproximación de Lanczos, precisa en todo el eje positivo y económica de evaluar. Debido a que cada cantidad grande se mantiene como su logaritmo hasta la 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 log-gamma único en lugar del factorial en el denominador. HotXLS devuelve 0.2460938 para BINOM.DIST(5,10,0.5,FALSE) and 0.6766764 para la acumulativa POISSON.DIST(2,2,TRUE), coincidiendo con Excel en todos los dígitos que imprime.

Inversas acotando la curva directa

Una función de distribución inversa hace la pregunta opuesta: dada una probabilidad, encontrar el valor cuya CDF sea igual a ella. Solo una inversa en este conjunto tiene 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 a la precisión de un double en todo el rango, dividida 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 elegido del soporte de la distribución, luego realiza una 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 gamma y chi-cuadrado, el acotamiento comienza en cero y 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 t acota límites simétricos que se ensanchan hacia afuera; la inversa F realiza una bisección en un intervalo no negativo. El costo es de unas pocas docenas de evaluaciones de CDF por llamada, lo cual es invisible a la velocidad de la hoja de cálculo, y el beneficio es que cada inversa es exactamente tan precisa como la función directa que invierte. Por eso, un viaje de ida y vuelta como CHISQ.DIST(CHISQ.INV(0.7,5),5,TRUE) devuelve 0.7 a una diferencia insignificante.

El logaritmo en base 10 que se ocultó en la cola

Este es el error que vale la pena contar, porque es de los que sobreviven mucho tiempo. La rutina normal inversa de Acklam tiene tres ramas. La rama central ancha, utilizada siempre que la probabilidad se sitúa entre aproximadamente 0.025 y 0.975, pasa la entrada a través de una relación polinómica sin ningún logaritmo en ella. Las dos ramas de la cola, para probabilidades muy pequeñas o muy grandes, toman primero un logaritmo de la entrada, porque la cola se comporta como la raíz cuadrada de menos el logaritmo natural de p.

Una versión temprana de la rama de la cola tomaba un logaritmo en base 10 donde correspondía el logaritmo natural. Los dos difieren por un factor constante de aproximadamente 2.30, por lo que los resultados de la cola eran incorrectos por un margen constante y considerable. And, sin embargo, la función se veía bien en cada verificación casual, porque las verificaciones casuales viven en el medio. NORM.S.INV(0.5) es cero, NORM.S.INV(0.975) es el valor de libro de texto 1.959964, y ambos pasan por el polinomio central que nunca llama a un logaritmo en absoluto. 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 volvía desviado por la relación entre el logaritmo natural y el de base 10. Cualquier función que dependa de la normal inversa en su cola, incluidas las herramientas 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 rota poco común. La solución fue un cambio de 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 acumulativa t de Student conlleva una sutileza que es fácil de entender al revés. Su valor proviene de la 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 significa que la conversión depende de qué lado de cero caiga 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. Devuelva el valor beta directamente y reportará el lado incorrecto de la distribución, desviado por todo el cuerpo de la curva para cualquier x distinto de cero. Las variantes de cola derecha y de dos colas se basan en la misma rama, razón por la cual T.DIST.2T(1,1) vuelve como 0.5 y T.DIST(1,1,TRUE) como 0.75, y la inversa T.INV realiza una bisección contra esta CDF corregida para que el viaje de ida y vuelta se complete.

Nada de esto es visible desde la celda, y ese es el resultado previsto. Escribe una fórmula y lee un número que coincide con Excel. Si está extendiendo el motor con su propia lógica, la mecánica de registrar una función se cubre en nuestro tutorial sobre el motor de fórmulas y funciones personalizadas, y la forma en que las fórmulas llegan a través de hojas y rangos con nombre 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, creación de gráficos y formato cubiertas en otras partes de este blog.