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

Escriba =NORM.DIST(115,100,15,TRUE) en una celda y Excel devuelve 0.8413447 sin ceremonia. La llamada parece una simple búsqueda. No lo es. Detrás de ese único número se encuentra la distribución normal acumulativa, una integral sin forma cerrada, y detrás de CHISQ.INV.RT y BETA.DIST se ocultan funciones especiales que una biblioteca cuidadosa tiene que evaluar, no aproximar a mano. Un componente de hoja de cálculo que afirma tener compatibilidad con Excel tiene que reproducir estos valores hasta el último dígito que Excel muestra, 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 totalmente 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 bifurcación que mantienen estable la aritmética y un error en la normal inversa que se ocultó en la cola durante mucho tiempo porque el caso común nunca tocaba la línea defectuosa

Una llamada de hoja de cálculo, cincuenta distribuciones detrás de ella

Las funciones abarcan las familias que un libro de estadística suele buscar. 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 los ayudantes de inferencia como CONFIDENCE.T y CONFIDENCE.NORM. Desde la posición del que llama, cada una es una única fórmula. Usted establece las entradas en las celdas, le pide al libro 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;   // observación
  sh.Range['A2', 'A2'].Value := 100;   // media
  sh.Range['A3', 'A3'].Value := 15;    // desviación estándar

  // El analizador de fórmulas XLS usa ';' como separador de argumentos.
  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 en el que la gente suele tropezar 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 celda almacenadas 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 estas funciona en Calculate, el resto sigue la misma ruta

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

La mayoría de las distribuciones acumulativas en este conjunto no se calculan sumando o integrando sus propias definiciones. Se calculan a partir de dos funciones especiales: la gamma incompleta inferior regularizada, escrita P(a, x), y la beta incompleta regularizada, escrita I_x(a, b). Internamente, estos son los ayudantes en los que se apoyan los despachadores, 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 beta incompleta regularizada en los argumentos correctos. La CDF de Poisson es la gamma incompleta superior Q. Implemente bien las funciones gamma y beta y una docena de distribuciones heredarán su precisión gratis

La palabra "regularizada" es el punto principal. La gamma incompleta cruda crece como un factorial y la integral beta cruda puede sufrir desbordamiento (underflow o overflow) mucho antes de que lo haga la respuesta. Las formas regularizadas se dividen por la gamma o beta completa, por lo que viven enteramente 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 a un chi-cuadrado con dos grados de libertad y a uno con doscientos sin que los términos intermedios se salgan del final 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 conlleva 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 de ella

La rutina gamma incompleta toma una decisión antes de calcular algo: compara x contra a + 1. Ese límite no es arbitrario. La expansión de la serie de potencias de P(a, x) converge rápidamente cuando x es pequeña en relación con a, y lentamente, eventualmente inútilmente, cuando x es grande. La fracción continua tiene el carácter opuesto. Por lo tanto, el motor usa 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 piso diminuto: cada vez que un término intermedio cae por debajo de aproximadamente 1e-30 en magnitud, se limita a 1e-30, lo que mantiene la recurrencia finita sin alterar el valor convergido. La misma limitación aparece en la fracción continua de la beta incompleta por la misma razón. Es una constante pequeña haciendo un trabajo 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 como máximo k eventos con media λ es Q(k + 1, λ). Enrutarlo 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 muchas 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 cincuenta y dos sobre veintiséis es un número entero enorme. Fórmelo directamente y el numerador desbordará un double antes de la división que lo devolvería a una probabilidad sensata. El motor nunca lo 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 exponencia una vez al final

// Masa de probabilidad binomial, evaluada completamente en espacio logarítmico.
// LnGammaF(n+1) es ln(n!); los tres log-factoriales forman ln(C(n,k)),
// y todo el exponente se construye antes de una sola llamada a Exp.
//   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 barata de evaluar. Dado que cada cantidad grande se mantiene como su logaritmo hasta el Exp final, el número más grande que la rutina materializa alguna vez es la probabilidad en sí, que es como máximo uno. La función de masa de Poisson sigue la misma receta, con el único término log-gamma sustituyendo al factorial en el denominador. Las formas cerradas son casos especiales en los bordes, donde p es exactamente cero o uno, por lo que el código nunca llama a Ln(0). HotXLS devuelve 0.2460938 para BINOM.DIST(5,10,0.5,FALSE) y 0.6766764 para el POISSON.DIST(2,2,TRUE) acumulativo, coincidiendo con Excel en todos los dígitos que imprime

Inversas acotando la curva hacia adelante

Una función de distribución inversa plantea la pregunta opuesta: dada una probabilidad, encuentre el valor cuya CDF la iguala. Solo una inversa en este conjunto tiene una fórmula directa rápida. NORM.S.INV, la normal estándar inversa, utiliza una aproximación racional de Acklam, un par de razones polinómicas precisas aproximadamente a la precisió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 del soporte de la distribución, luego biseca: evalúa la CDF hacia adelante en el punto medio, mueve cualquier límite que mantenga encerrada la probabilidad objetivo y repite hasta que el intervalo sea estrecho. Para las inversas de gamma y chi-cuadrado, el corchete 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 biseca en un intervalo no negativo. El costo es de unas pocas docenas de evaluaciones de CDF por llamada, que 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 hacia adelante que invierte. Es por eso que un viaje de ida y vuelta como CHISQ.DIST(CHISQ.INV(0.7,5),5,TRUE) devuelve 0.7 por un pelo

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

Aquí está el error que vale la pena contar, porque es del tipo que sobrevive mucho tiempo. La rutina inversa-normal de Acklam tiene tres ramas. La rama central ancha, que se usa siempre que la probabilidad se encuentra entre aproximadamente 0.025 y 0.975, pasa la entrada a través de una razón polinómica sin logaritmo en ninguna parte. 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 pertenecía el logaritmo natural. Los dos difieren por un factor constante de aproximadamente 2.30, por lo que los resultados de la cola estaban equivocados por un margen consistente y considerable. Y 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 1.959964 de libro de texto, y ambos corren a través del polinomio central que nunca llama a un logaritmo en absoluto. El error solo apareció una vez que una probabilidad cruzó hacia una cola, digamos NORM.S.INV(0.001), que debe devolver -3.0902323 y en su lugar regresaba con la diferencia de proporción de log-natural-versus-base-10. Cualquier función que dependa de la normal inversa en su cola, incluidos los ayudantes de intervalo de confianza, heredó el mismo sesgo. La lección es mundana y costosa: una función con una estructura de rama necesita puntos de prueba dentro de cada rama, porque un camino común correcto enmascarará alegremente uno raro roto. La corrección fue un cambio de un solo token del logaritmo 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 tiene una sutileza de la que es fácil equivocarse. 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 acumulativa hasta x. La forma simétrica de la distribución t significa que la conversión depende del lado de cero en el que caiga x

// CDF de t de Student. ib es la beta incompleta regularizada en df/(df+x*x),
// que mide la cola simétrica. El valor acumulativo depende de
// el signo de x; devolver ib sin convertir da la cola incorrecta.
ib := BetaIF(df / 2, 0.5, df / (df + x * x));
if x > 0 then
  result := 1 - 0.5 * ib        // por encima de la media: uno menos la mitad de la cola
else if x < 0 then
  result := 0.5 * ib            // por debajo de la media: la mitad de la cola
else
  result := 0.5;                // exactamente en la media

Para x por encima de cero, la probabilidad acumulativa 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 equivocado de la distribución, desfasado por todo el cuerpo de la curva para cualquier x distinto de cero. Las variantes de cola derecha y dos colas se basan en la misma rama, por lo que T.DIST.2T(1,1) vuelve como 0.5 y T.DIST(1,1,TRUE) como 0.75, y la inversa T.INV biseca contra esta CDF corregida para que el viaje de ida y vuelta se cierre

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 concuerda con Excel. Si está ampliando el motor con su propia lógica, la mecánica de registro de una función se trata en nuestro recorrido por el motor de fórmulas y las funciones personalizadas, y la forma en que las fórmulas llegan a otras hojas y rangos con nombre se trata en el artículo sobre nombres definidos y fórmulas entre hojas. Todo esto se envía 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 que se tratan en otros lugares de este blog