¿Por qué es tan importante la familia exponencial en estadística?

Question

¿Por qué es tan importante la familia exponencial en estadística?

Hace poco leí sobre la familia exponencial dentro de las estadísticas. Según tengo entendido, la familia exponencial se refiere a cualquier función de distribución de probabilidad que se puede escribir en el siguiente formato (nótese el "exponente" en esta ecuación):

Esto incluye funciones de distribución de probabilidad comunes, como la distribución normal El distribución gamma El Distribución de Poisson etc. Las distribuciones de probabilidad de la familia exponencial se utilizan a menudo como "función de enlace" en los problemas de regresión (por ejemplo, en entornos de datos de recuento, la variable de respuesta puede relacionarse con las covariables a través de una distribución de Poisson) - las funciones de distribución de probabilidad que pertenecen a la familia exponencial se utilizan a menudo debido a sus "propiedades matemáticas deseables". Por ejemplo, estas propiedades son las siguientes:

¿Por qué son tan importantes estas propiedades?

A) La primera propiedad se refiere a la "estadística suficiente". Una "estadística suficiente" es una estadística que proporciona más información para cualquier conjunto de datos/parámetro del modelo en comparación con cualquier otra estadística.

Me cuesta entender por qué esto es importante. En el caso de la regresión logística, se utiliza la función de enlace logit (parte de la familia exponencial) para vincular la variable de respuesta con las covariables observadas. ¿Qué son exactamente los "estadísticos" en este caso (por ejemplo, en un modelo de regresión logística, estos "estadísticos" se refieren a la "media" y la "varianza" de los coeficientes beta del modelo de regresión)? ¿Cuáles son los "valores fijos" en este caso?

B) Las familias exponenciales tienen priores conjugados.

En el entorno bayesiano, una prioridad p(thetha | x) se llama prioridad conjugada si está en la misma familia que la distribución posterior p(x | thetha). Si una a priori es una a priori conjugada, esto significa que existe una solución de forma cerrada y que las técnicas de integración numérica (por ejemplo, MCMC ) no son necesarios para muestrear la distribución posterior. ¿Es esto correcto?

C) ¿La tercera propiedad es esencialmente similar a la segunda?

D) No entiendo en absoluto la cuarta propiedad. Bayas Variacionales son una alternativa a las técnicas de muestreo MCMC que aproximan la distribución posterior con una distribución más simple - esto puede ahorrar tiempo computacional para distribuciones posteriores de alta dimensión con grandes datos. ¿Significa la cuarta propiedad que los Bayes variacionales con priores conjugados en la familia exponencial tienen soluciones de forma cerrada? Así que cualquier modelo bayesiano que utilice la familia exponencial no requiere MCMC - ¿es esto correcto?

Referencias:

• _Familia exponencial_
• _Estadística suficiente_
• _Conjugado previo_

Answer 1

Excelentes preguntas.

Con respecto a A: Una estadística suficiente no es más que una destilación de la información que contiene la muestra con respecto a un modelo determinado . Como es de esperar, si tiene una muestra $x_i \sim N(\mu,\sigma^2)$ para $i \in \{1, \ldots, N\}$ y cada una independiente, está claro que mientras calculemos la media y la varianza muestral, no importa cuáles sean los valores de cada $x_i$ son. En la regresión lineal (más fácil de hablar que en la logística en este contexto), la distribución muestral del vector de coeficientes desconocidos (para una varianza conocida) es $N(\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\mathbf{y}, \sigma^2\mathbf{X}^\top\mathbf{X})^{-1})$ por lo que mientras estas cantidades finales sean idénticas, la inferencia basada en ellas también lo será. Esta es la idea de suficiencia.

Tenga en cuenta que en el $N(\mu,\sigma^2)$ ejemplo, la estadística suficiente se compone de sólo dos números: $\hat{\mu}=\frac{1}{N}\sum_{i=1}^N x_i$ y $\frac{1}{N}\sum_{i=1}^N (x_i-\hat{\mu})^2$ No importa cuán grande sea el tamaño de nuestra muestra $N$ es (y suponiendo que $N>2$ ). Asimismo, el vector $(\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\mathbf{y}$ es de dimensión $P$ y $\sigma^2(\mathbf{X}^\top\mathbf{X})^{-1}$ de dimensión $P\times P$ (aquí $P$ es la dimensión de la matriz de diseño), que son independientes de $N$ (aunque, técnicamente, la matriz $\sigma^2(\mathbf{X}^\top\mathbf{X})^{-1}$ es sólo una constante bajo nuestros supuestos). Así que en estos ejemplos, la estadística suficiente tiene un número de valores (no valores fijos ), o como yo diría, de dimensión fija.

Anotemos tres cosas más. En primer lugar, que no existe tal cosa como el estadística suficiente para una distribución, sino que hay muchas estadísticas posibles que pueden ser suficientes, y que pueden ser de diferente dimensión. De hecho, lo segundo que tenemos que discutir es que la muestra completa en sí misma, ya que naturalmente contiene toda la información que contiene, es siempre una estadística suficiente. Se trata de un caso trivial, pero importante, ya que en general no se puede esperar encontrar siempre una estadística suficiente de dimensión inferior a $N$ . Y la última cosa a tener en cuenta es la especificidad del modelo: por eso escribí con respecto a un determinado modelo anterior . El cambio de la probabilidad cambiará las estadísticas suficientes, al menos potencialmente, para un conjunto de datos determinado.

Con respecto a B : Lo que dices es correcto, pero además de permitir posteriors analíticos en el caso univariante, la conjugación tiene serias ventajas en el contexto de los modelos jerárquicos bayesianos estimados mediante MCMC. Esto se debe a que condicional también están disponibles en forma cerrada. Así que podemos acelerar los algoritmos MCMC del estilo Metrópolis dentro de Gibbs con la conjugación.

Con respecto a C: Es definitivamente una idea similar, pero quiero dejar claro que estamos hablando de dos distribuciones diferentes: "posterior" versus "posterior predictiva". Como su nombre indica, ambas son distribuciones posteriores, lo que significa que son distribuciones de una variable desconocida condicionadas a nuestros datos conocidos. Una "posterior" simple y llanamente se refiere a algo como $P(\mu, \sigma^2| \{x_1, \ldots, x_N\})$ de nuestro ejemplo normal anterior: una distribución de parámetros desconocidos definida en las distribuciones generadoras de datos. En cambio, una "predicción posterior" da la distribución de un hipotético $N+1$ 'st punto de datos $x_{N+1}$ condicionada a los datos observados: $P(x_{N+1}| \{x_1, \ldots, x_N\})$ . Observe que esto no está condicionado a los parámetros $\mu$ y $\sigma^2$ : había que integrarlos fuera. Es esta integral adicional la que está garantizada por la conjugación.

Con respecto a D: En el contexto de Bayes Variacional (VB), se tiene una distribución posterior $P(\theta|X)$ donde $\theta$ es un vector de $P$ parámetros y $X$ son algunos datos. En lugar de intentar generar una muestra a partir de ellos, como hace MCMC, vamos a utilizar una distribución posterior aproximada con la que es fácil trabajar y que se acerca bastante a la verdadera. Eso se llama una distribución variacional y se denota $Q_\eta(\theta)$ . Obsérvese que nuestra distribución variacional está indexada por parámetros variacionales $\eta$ . Los parámetros variacionales no se parecen en nada a los parámetros sobre los que hacemos la inferencia bayesiana, y no se parecen en nada a nuestros datos. No tienen una distribución asociada y no tienen un papel hipotético generando los datos. Más bien, se eligen como resultado de un algoritmo de optimización iterativo. La idea de la inferencia variacional es definir alguna medida de disimilitud entre la distribución variacional y la verdadera posterior y luego minimizar esa medida con respecto a los parámetros $\eta$ . Denotaremos el resultado de ese proceso de optimización por $\hat{\eta}(X)$ . En ese momento, con suerte $Q_{\hat{\eta}(X)}(\theta)$ está bastante cerca de $P(\theta|X)$ y si hacemos inferencias utilizando $Q_{\hat{\eta}(X)}(\theta)$ en su lugar obtendremos respuestas similares.

Ahora bien, ¿dónde encaja la conjugación? Una forma popular de medir la disimilitud es esta medida, que se llama KL inverso costo:

$$\hat{\eta}(X) := \underset{\eta}{\textrm{argmin}}\, \mathbb{E}_{\theta\sim Q_\eta}\bigg[\frac{\log Q_{\eta}(\theta)}{\log P(\theta|X)}\bigg]$$

Esta integral no puede resolverse en términos de funciones simples en general. Sin embargo, está disponible en forma cerrada cuando:

1. Utilizamos un previo conjugado para definir $P(\theta|X)$ .

2. Suponemos que la distribución variacional es independiente En otras palabras, que $q_\eta(\theta)=\prod_{j=1}^P q_{j,\eta}(\theta_j)$ .

3. Además, nos limitamos a una $q_{j,\eta_j}$ para cada $j$ (que se determina por la probabilidad).

Así que no es que la variante posterior esté disponible en forma cerrada. Más bien, es que la función de coste que define la variante posterior está disponible en forma cerrada. El hecho de que la función de coste sea de forma cerrada hace que el cálculo de la distribución variacional sea un problema de optimización más fácil, ya que podemos calcular analíticamente los valores de la función y los gradientes.

Answer 2

Dado que no veo el resultado mencionado en ese hilo hasta ahora, permítanme mencionar una cuestión a menudo descuidada aunque significativa, a saber, que, en relación con la pregunta A Las familias exponenciales están estrechamente relacionadas con la noción de suficiencia debido a la Lema de Pitman-Koopman-Darmois :

Supongamos que $X_n$ , $n = 1 , 2 , 3 , \dots$ son idénticos e independientes distribuidas aleatorias cuya distribución se sabe que está en alguna familia de distribuciones de probabilidad con soporte fijo. Sólo si esa familia es una familia exponencial existe un estadístico suficiente (posiblemente de valor vectorial) $\displaystyle T(X_{1},\dots ,X_{n})$ cuyo número de componentes escalares no aumenta con el tamaño de la muestra $n$ aumenta. Wikipedia

En resumen, no hay ninguna familia de soporte fijo fuera de las familias exponenciales que goce de un estadístico suficiente de dimensión fija.

Y ya que estoy en ello, apunta B no es un argumento convincente, ya que las expresiones de forma cerrada no son un argumento sólido para definir una distribución a priori que refleje las creencias previas de uno. Además, existen familias exponenciales cuyos conjugados no son manejables, como la distribución Beta. Por lo tanto, los conjugados se ven mejor como aproximaciones útiles en el espíritu de la pregunta D .

Answer 3

Grandes preguntas. Hay muchas maneras de responderlas. John Madden hace un excelente trabajo, pero voy a citar un poco la respuesta de Ben aquí con respecto a las estadísticas suficientes.

La función de pérdida para un modelo lineal gaussiano (como señala Ben) es

$$\ell_{\mathbf{y}, \mathbf{x}}(\boldsymbol{\beta}, \sigma)=-n \ln \sigma-\frac{1}{2 \sigma^{2}}\|\mathbf{y}-\mathbf{x} \boldsymbol{\beta}\|^{2}$$

En mucho material pedagógico sobre el ajuste de modelos, utilizaríamos los datos para calcular los gradientes de la pérdida y realizar algún tipo de rutina de optimización. Un código para hacer esto podría ser como

def compute_loss(beta, X, y):

    number_rows = len(X)
    loss = 0

    # Potentially expensive!
    for i in range(number_rows):
        loss += (y[i] - X[i]@beta)**2

    loss/=number_rows

    return loss

Si tenemos muchos datos, entonces este bucle (y cualquier iteración adicional sobre los datos para calcular los gradientes, por ejemplo) podría ser costoso de calcular. Pero, como podemos escribir suficientes estadísticas para la familia exponencial, podemos mejorar nuestro cálculo drásticamente.

Como escribe Ben, la pérdida se puede reescribir como

$$\begin{aligned} &=-n \ln \sigma-\frac{1}{2 \sigma^{2}} \mathbf{y}^{\mathrm{T}} \mathbf{y}-\frac{1}{2 \sigma^{2}}\left(2 \boldsymbol{\beta}^{\mathrm{T}} \mathbf{T}_{1}-\boldsymbol{\beta}^{\mathrm{T}} \mathbf{T}_{2} \boldsymbol{\beta}\right) \\ &=h(\mathbf{y}, \sigma)+g_{\beta}\left(\mathbf{T}_{1}, \mathbf{T}_{2}, \sigma\right) \end{aligned}$$

Donde $\mathbf{T}_{1} \equiv \mathbf{T}_{1}(\mathbf{x}, \mathbf{y}) \equiv \mathbf{x}^{\mathrm{T}} \mathbf{y}$ y $\mathbf{T}_{2} \equiv \mathbf{T}_{2}(\mathbf{x}, \mathbf{y}) \equiv \mathbf{x}^{\mathrm{T}} \mathbf{x}$ . Así, podemos calcular estas cantidades una sola vez y ajustar nuestros modelos, en lugar de volver a calcularlas en cada paso de actualización.

Comparemos en un problema simulado. Asumiré $\sigma=1$ y no requiere estimación por simplicidad. Variaré el número de observaciones y el número de covariables en el modelo y compararé el tiempo para optimizar la pérdida asumiendo cuando se utilizan estadísticas suficientes y cuando se utiliza un enfoque ingenuo. A continuación se muestra un gráfico del tiempo esperado hasta la finalización más/menos una desviación estándar. Podemos ver que a medida que los datos se hacen más grandes en tamaño, el uso de estadísticas suficientes es ventajoso.

Ahora bien, toda esta respuesta es un poco un hombre de paja. No te he mostrado los resultados cuando no se usan estadísticas suficientes y no se usan bucles (así que quizás aprovechando algo de álgebra lineal). Pero el punto sobre la estadística suficiente se mantiene. Que podamos representar toda la información de una muestra con un solo número, calculado una vez, es una propiedad muy valiosa.

Código para reproducir los experimentos:

import numpy as np
import pandas as pd
from scipy.optimize import minimize
from itertools import product
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_theme(style="darkgrid")

# First, simulate big data!

def make_data(N, p):
    X = np.random.normal(size = (N,p))
    beta = np.random.normal(2, 2, size = p)
    y = X@beta + np.random.normal(size = N)
    return X,y

# Next, set up a loss function to optimize using loops
def solve_naive(N, p):
    X, y = make_data(N, p)
    def loss_and_grad(w):

        number_rows, number_columns = X.shape
        grads = np.zeros_like(w)
        loss=0
        for i in range(number_rows):
            res= (y[i] - X[i]@w)
            loss+= res**2/number_rows

            grads+= -2*res*X[i]/number_rows

        return loss, grads

    r=minimize(loss_and_grad, x0=np.zeros(p), jac=True)
    return loss_and_grad

# Next, set up a loss function to opimize which only uses sufficient statistics
def solve_sufficient(N, p):
    X, y = make_data(N, p)

    T1 = X.T@y
    T2 = X.T@X
    const = y.T@y
    n = len(X)

    def loss_and_grad(w):

        loss = const - (2*w@T1 - w@T2@w)

        grads = -2*T1 + (T2@w + T2.T@w)

        return loss/n, grads/n

    r=minimize(loss_and_grad, x0=np.zeros(p), jac=True)
    return loss_and_grad

# A helper function to time the optimization for various datasets
def time_to_optimize(N, p):

    naive_optimization_times = %timeit -o -n 10 solve_naive(N, p)
    suff_optimization_times = %timeit -o -n 10 solve_sufficient(N, p)

    suff = pd.DataFrame({'times':suff_optimization_times.timings,
                         'type': 'Sufficient Statistics'})

    naive = pd.DataFrame({'times':naive_optimization_times.timings,
                         'type': 'Naive'})

    df = pd.concat((suff, naive))

    df['N'] = N
    df['p'] = p

    return df

if __name__ == '__main__':

    Ns = [1_000, 10_000, 100_000]
    ps = [10, 100, 250]

    prods = product(Ns, ps)
    frames = []

    for N, p in prods:

        frames.append(time_to_optimize(N, p))

    df = pd.concat(frames).reset_index(drop=True)

    fig, ax = plt.subplots(dpi = 240, figsize = (8, 5))
    grid = sns.lineplot(data=df, x='N', y='times', hue='type', style='p')

    grid.set(xscale="log", yscale="log", xlabel = 'Number of Observations', ylabel='Execution Time (Seconds)')
    grid.legend(loc='best', prop = {'size':6})

Answer 4

Ninguna de las propiedades que menciona OP era importante cuando el más popular Las distribuciones de la familia exponencial fueron descubiertas o puestas en uso. Es decir no decir que las propiedades son irrelevantes o no son importantes. Todas estas características son interesantes y útiles, pero son no las razones de la popularidad de las distribuciones. Distribuciones como la Gaussiana y la Poisson se hicieron populares porque son omnipresentes en la naturaleza, y también se aproximan muy bien a los procesos en la ingeniería y la industria.

Poisson. He aquí un ejemplo: La distribución de Poisson. Ésta y sus distribuciones afines, como la exponencial, se observan en el fenómeno de la desintegración radiactiva. No sólo es simplemente conveniente aproximar los datos observados con ellas, sino que surgen de la teoría de la desintegración radiactiva en su forma exacta, y las observaciones son consistentes con las ecuaciones más allá de cualquier duda con una precisión arbitraria limitada sólo por la precisión de nuestra instrumentación. Esto ha sido descubierto a principios del siglo XX y la precisión cada vez mayor de las mediciones no ha producido ninguna desviación de las ecuaciones exactas. En este punto es simplemente un hecho que la distribución de Poisson explica el número de núcleos que decaen en un periodo de tiempo determinado, y que el tiempo entre dos eventos es de distribución exponencial.

Gaussiano. La distribución gaussiana se observa en numerosos fenómenos naturales, como el movimiento browniano. He puesto un ejemplo exótico en mi comentario a la pregunta: El principio de incertidumbre de Heisenberg (UP) de la mecánica cuántica (QM).

Lo primero que hay que señalar aquí es que los físicos no han encontrado hasta ahora ninguna desviación de la conclusión de la QM. La precisión con la que se ha confirmado que funciona es alucinante y no tiene precedentes en ningún campo de estudio de las ciencias. Por lo tanto, se puede aceptar con seguridad la UP como un hecho de la naturaleza.

En su forma simple establece que las incertidumbres sobre la medición del momento y la coordenada de una partícula simultáneamente tiene un límite inferior: $$\Delta p\Delta x\ge h/4\pi$$ Si se mide la coordenada con mucha precisión, la medición del momento será muy imprecisa, etc.

Pero cuando la desigualdad se convierte en igualdad $\Delta p\Delta x= h/4\pi$ ? Para un paquete de ondas gaussianas.

No siempre tenemos elección. En los ejemplos anteriores no elegimos Poisson o Gaussian porque tienen algunas buenas propiedades estadísticas. Los utilizamos porque debe utilizarlos. Así es como funciona la naturaleza, y si queremos describir los fenómenos observados, no tenemos más remedio que utilizar estas distribuciones, independientemente de su conveniencia o elegancia.