Más sobre ciencia de datos y estadística en: cienciadedatos.net

Introducción

Identificar el tipo de distribución teórica que mejor se ajusta a una variable empírica es un paso crítico en el análisis estadístico. Es la base para realizar contrastes de hipótesis paramétricos, estimar intervalos de confianza y preprocesar datos en modelos de aprendizaje estadístico y Machine Learning (especialmente en modelos lineales o gaussianos).

En Python, el ecosistema de análisis de datos ofrece herramientas potentes para esta tarea. Este documento explora las funcionalidades del módulo scipy.stats, enfocándose en cómo ajustar y comparar múltiples distribuciones para determinar cuál describe mejor los datos observados.

Ajuste de una distribución

Ajustar una distribución paramétrica consiste en determinar los valores de sus parámetros de forma que la distribución teórica represente lo mejor posible el comportamiento de un conjunto de datos reales. Una vez estimados estos parámetros, la distribución queda completamente definida y puede usarse para describir, analizar o predecir el comportamiento del sistema estudiado. Por ejemplo, una distribución normal queda totalmente definida por dos parámetros: su ubicación (media, $\mu$) y su escala (desviación estándar, $\sigma$).

Existen varios métodos que permiten encontrar los parámetros óptimos, uno de los más utilizados, y el que implementa scipy.stats, es el método de Estimación de Máxima Verosimilitud (MLE, por sus siglas en inglés: Maximum Likelihood Estimation). Este método busca los parámetros que maximizan la probabilidad conjunta de observar los datos disponibles.

La librería scipy.stats dispone de más de 90 distribuciones implementadas. Puedes consultar el listado completo en su documentación oficial:

La definición anterior pretende ser intuitiva. Formalmente, el ajuste de una distribución paramétrica se puede definir de la siguiente manera:

El ajuste de una distribución paramétrica consiste en seleccionar una familia paramétrica de distribuciones $\{f(x;\theta), \theta \in \Theta\}$ y estimar el valor del parámetro desconocido $\theta$ a partir de una muestra aleatoria $X_1, \dots, X_n$ de una variable aleatoria $X$, bajo el supuesto de que la distribución de $X$ pertenece a dicha familia. La estimación de $\theta$ se realiza mediante un criterio estadístico formal, como el método de máxima verosimilitud, el método de los momentos o el de mínimos cuadrados, de modo que la distribución $f(x;\hat{\theta})$ resulte óptima con respecto a dicho criterio.

Además de diferenciar entre distribuciones continuas y discretas, es útil poder seleccionarlas por el soporte, es decir, el conjunto de valores sobre los que la distribución está definida. Por ejemplo, si se quiere modelar la velocidad del viento, aunque no se conozca el tipo exacto de distribución, se puede restringir la selección a aquellas distribuciones cuyo soporte esté contenido en $[0,+\infty)$.

# Distribuciones disponibles en scipy.stats
# ==============================================================================
from scipy import stats
import pandas as pd
import numpy as np

distribuciones = [
    getattr(stats, d) for d in dir(stats)
    if isinstance(getattr(stats, d), (stats.rv_continuous, stats.rv_discrete))
]

data = []
for dist in distribuciones:
    try:
        a = getattr(dist, 'a', np.nan)
        b = getattr(dist, 'b', np.nan)
        name = getattr(dist, 'name', str(dist))
        data.append({'distribucion': name, 'soporte_min': a, 'soporte_max': b})
    except Exception as e:
        continue

info_distribuciones = pd.DataFrame(data)
info_distribuciones = info_distribuciones.sort_values(by=['soporte_min', 'soporte_max']).reset_index(drop=True)
display(info_distribuciones)
distribucion soporte_min soporte_max
0 levy_l -inf 0.0
1 weibull_max -inf 0.0
2 cauchy -inf inf
3 crystalball -inf inf
4 dgamma -inf inf
... ... ... ...
127 pareto 1.0 inf
128 truncpareto 1.0 inf
129 yulesimon 1.0 inf
130 zipf 1.0 inf
131 zipfian 1.0 inf

132 rows × 3 columns

Comparación de ajustes

El método de ajuste por Maximum Likelihood Estimation (MLE) permite, dada una distribución específica, encontrar los valores de sus parámetros que mejor se ajustan a los datos.

En la práctica, a menudo no se conoce de antemano la distribución que siguen los datos. Con experiencia, el analista puede acotar las distribuciones candidatas (continuas o discretas, solo positivas, etc.), pero es necesario contar con un criterio que cuantifique la bondad de ajuste de cada distribución y permita comparar entre ellas. Un libro en el que se explican las características de las principales distribuciones probabilísticas es Distributions for Modeling Location, Scale, and Shape Using GAMLSS in R By Robert A. Rigby, Mikis D. Stasinopoulos, Gillian Z. Heller, Fernanda De Bastiani.

Generalized Akaike information criterion (GAIC)

En distribuciones paramétricas que han sido ajustadas siguiendo una estrategia de Maximum Likelihood Estimation (MLE), una forma de cuantificar la bondad de ajuste del modelo, es decir, lo bien que se ajusta a los datos, es empleando el propio valor de likelihood conseguido en el ajuste. ¿Qué significa esto?

El valor de likelihood devuelto por una distribución para una observación $x_i$ es una medida de la probabilidad con la que esa distribución podría haber generado dicha observación. Para facilitar los cálculos matemáticos, normalmente se utiliza el logaritmo ( log likelihood), pero la interpretación es la misma, a mayor log likelihood mayor probabilidad.

Si se suma el log likelihood de todas las observaciones con las que se ha ajustado la distribución, se dispone de un valor que mide la "compatibilidad" entre la distribución y los datos. Siguiendo esta idea es como se define el término fitted global deviance (GDEV), que no es más que el log likelihood del modelo multiplicado por $-2$.

$$\text{GDEV} = − 2 \log(likelihood)$$

El problema de emplear el GDEV para comparar distribuciones es que no tiene en cuenta los grados de libertad de cada una (su flexibilidad). En términos generales, cuantos más parámetros tenga una distribución, con más facilidad se ajusta a los datos y mayor es su log likelihood. Esto significa que utilizar el GDEV o el log likelihood para comparar distribuciones con distinto número de parámetros no es una estrategia justa, casi siempre ganará la que más parámetros tenga aunque no sea realmente la que mejor describe el comportamiento de los datos.

Una forma de mitigar este problema es mediante el uso del GAIC (generalized Akaike information criterion), que incorpora una penalización $k$ por cada parámetro que tenga la distribución:

$$\text{GAIC} = \text{GDEV} + k \times \text{nº parámetros}= $$$$= −2 \log(likelihood) + k \times \text{nº parámetros}$$

Dependiendo del grado de penalización, se favorece más o menos la sencillez del modelo. Dos estándares de esta métrica son el AIC (Criterio de información de Akaike) y BIC (Bayesian information criterion) también conocida como SBC. El primero utiliza como valor de penalización $k=2$ y el segundo $k=\log(\text{nº observaciones)}$.

$$\text{AIC} = −2 \log(\text{likelihood}) + 2 \times \text{nº parámetros}$$$$\text{BIC} = −2 \log(\text{likelihood}) + log (\text{nº observaciones}) \times \text{nº parámetros}$$

En la práctica, AIC suele favorecer distribuciones con un más parámetros (overfitting), mientras BIC/SBC tiende al contrario underfitting. Los autores del paquete del libro Distributions for Modeling Location, Scale, and Shape Using GAMLSS in R recomiendan el uso de valores de $k$ en el rango $2.5 \leq k \leq 4$ o de $k = \sqrt{\log(n)}$ cuando $n \geq 1000$.

La multiplicación por -2 convierte la maximización del log likelihood en un problema de minimización, de modo que valores más bajos de GDEV, AIC o BIC indican mejor ajuste.

Es importante tener en cuenta que, ninguna de estas métricas, sirven para cuantificar la calidad del ajuste en un sentido absoluto, sino para comparar la calidad relativa entre distribuciones. Si todas las distribuciones ajustadas son malas, no proporcionan ningún aviso de ello.

Ejemplo

En este ejemplo se procede a ajustar dos distribuciones, normal y gamma, con el objetivo de modelizar la distribución del precio de venta de diamantes. Además de realizar los ajustes, se representan gráficamente los resultados y se calculan las métricas de bondad de ajuste AIC, BIC y Log-Likelihood con el objetivo de comparar e identificar el mejor que distribución se ajusta mejor.

Librerías

Las librerías utilizadas en este documento son:

# Tratamiento de datos
# ==============================================================================
import pandas as pd
import numpy as np
import seaborn as sns

# Gráficos
# ==============================================================================
import matplotlib.pyplot as plt
plt.style.use('fivethirtyeight')
plt.rcParams['lines.linewidth'] = 1.5
plt.rcParams['font.size'] = 8

# Ajuste de distribuciones
# ==============================================================================
from tqdm.auto import tqdm
from scipy import stats
import inspect
from statsmodels.distributions.empirical_distribution import ECDF

# Configuración warnings
# ==============================================================================
import warnings
from contextlib import nullcontext
warnings.filterwarnings('once')

Datos

Para esta demostración se emplean como datos el precio de los diamantes disponible en data set diamonds de la librería seaborn, en concreto, la columna price.

# Datos
# ==============================================================================
datos = sns.load_dataset('diamonds')
datos = datos.loc[datos.cut == 'Fair', 'price']

Ajuste de distribuciones

Dos de los primeros pasos a la hora de analizar una variable son: calcular los principales estadísticos descriptivos y representar las distribuciones observadas (empíricas). Si los datos se almacenan en un Serie de Pandas, pueden obtenerse los principales estadísticos descriptivos con el método describe().

# Estadísticos descriptivos
# ==============================================================================
datos.describe()

count     1610.000000
mean      4358.757764
std       3560.386612
min        337.000000
25%       2050.250000
50%       3282.000000
75%       5205.500000
max      18574.000000
Name: price, dtype: float64
# Gráficos distribución observada (empírica)
# ==============================================================================
fig, axs = plt.subplots(nrows=1, ncols=2, figsize=(8, 3))

# Histograma
axs[0].hist(x=datos, bins=30, color="#3182bd", alpha=0.5)
axs[0].plot(datos, np.full_like(datos, -0.01), '|k', markeredgewidth=1)
axs[0].set_title('Distribución empírica del precio diamantes')
axs[0].set_xlabel('precio')
axs[0].set_ylabel('counts')

# Función de Distribución Acumulada
# ecdf (empirical cumulative distribution function)
ecdf = ECDF(x=datos)
axs[1].plot(ecdf.x, ecdf.y, color="#3182bd")
axs[1].set_title('Función de distribución empírica')
axs[1].set_xlabel('precio')
axs[1].set_ylabel('CDF')

plt.tight_layout();

# Ajuste distribución normal
#===============================================================================
# 1) Se define el tipo de distribución
distribucion = stats.norm

# 2) Con el método fit() se obtienen los parámetros
parametros = distribucion.fit(data=datos)

# 3) Se crea un diccionario que incluya el nombre de cada parámetro
nombre_parametros = [p for p in inspect.signature(distribucion._pdf).parameters \
                     if not p=='x'] + ["loc","scale"]
parametros_dict = dict(zip(nombre_parametros, parametros))

# 3) Se calcula el log likelihood
log_likelihood = distribucion.logpdf(datos.to_numpy(), *parametros).sum()

# 4) Se calcula el AIC y el BIC
aic = -2 * log_likelihood + 2 * len(parametros)
bic = -2 * log_likelihood + np.log(datos.shape[0]) * len(parametros)

# 5) Gráfico
x_hat = np.linspace(min(datos), max(datos), num=100)
y_hat = distribucion.pdf(x_hat, *parametros)

fig, ax = plt.subplots(figsize=(6, 3))
ax.plot(x_hat, y_hat, linewidth=2, label=distribucion.name)
ax.hist(x=datos, density=True, bins=30, color="#3182bd", alpha=0.5)
ax.plot(datos, np.full_like(datos, -0.01), '|k', markeredgewidth=1)
ax.set_title('Distribución precio diamantes')
ax.set_xlabel('precio')
ax.set_ylabel('Densidad de probabilidad')
ax.legend();

#6) Información del ajuste
print('---------------------')
print('Resultados del ajuste')
print('---------------------')
print(f"Distribución   : {distribucion.name}")
print(f"Soporte        : {[distribucion.a, distribucion.b]}")
print(f"Parámetros     : {parametros_dict}")
print(f"Log likelihood : {log_likelihood}")
print(f"AIC            : {aic}")
print(f"BIC            : {bic}")

---------------------
Resultados del ajuste
---------------------
Distribución   : norm
Soporte        : [-inf, inf]
Parámetros     : {'loc': np.float64(4358.757763975155), 'scale': np.float64(3559.2807303891086)}
Log likelihood : -15449.966194325283
AIC            : 30903.932388650566
BIC            : 30914.700367566522

Se repite el proceso, pero esta vez con la distribución gamma.

# Ajuste distribución gamma
#===============================================================================
# 1) Se define el tipo de distribución
distribucion = stats.gamma

# 2) Con el método fit() se obtienen los parámetros
parametros = distribucion.fit(data=datos)

# 3) Se crea un diccionario que incluya el nombre de cada parámetro
nombre_parametros = [p for p in inspect.signature(distribucion._pdf).parameters \
                     if not p=='x'] + ["loc","scale"]
parametros_dict = dict(zip(nombre_parametros, parametros))

# 3) Se calcula el log likelihood
log_likelihood = distribucion.logpdf(datos.to_numpy(), *parametros).sum()

# 4) Se calcula el AIC y el BIC
aic = -2 * log_likelihood + 2 * len(parametros)
bic = -2 * log_likelihood + np.log(datos.shape[0]) * len(parametros)

# 5) Gráfico
x_hat = np.linspace(min(datos), max(datos), num=100)
y_hat = distribucion.pdf(x_hat, *parametros)

fig, ax = plt.subplots(figsize=(6, 3))
ax.plot(x_hat, y_hat, linewidth=2, label=distribucion.name)
ax.hist(x=datos, density=True, bins=30, color="#3182bd", alpha=0.5)
ax.plot(datos, np.full_like(datos, -0.01), '|k', markeredgewidth=1)
ax.set_title('Distribución precio diamantes')
ax.set_xlabel('precio')
ax.set_ylabel('Densidad de probabilidad')
ax.legend();

#6) Información del ajuste
print('---------------------')
print('Resultados del ajuste')
print('---------------------')
print(f"Distribución   : {distribucion.name}")
print(f"Soporte        : {[distribucion.a, distribucion.b]}")
print(f"Parámetros     : {parametros_dict}")
print(f"Log likelihood : {log_likelihood}")
print(f"AIC            : {aic}")
print(f"BIC            : {bic}")

---------------------
Resultados del ajuste
---------------------
Distribución   : gamma
Soporte        : [0.0, inf]
Parámetros     : {'a': np.float64(14.399103124683087), 'loc': np.float64(-6244.786483560789), 'scale': np.float64(726.3476320693849)}
Log likelihood : -15235.639828174484
AIC            : 30477.27965634897
BIC            : 30493.431624722904

Tanto la métrica AIC como la BIC coinciden en que la distribución gamma se ajusta mejor a los datos (valores de AIC y BIC más bajos). Esto se puede corroborar fácilmente con la inspección gráfica de los resultados.

En este caso, dado que los valores de precio solo pueden ser positivos y se tienen una notable cola derecha, la distribución gamma era una candidata mucho mejor que la normal. Sin embargo, hay otras posibles distribuciones, algunas de las cuales podrían ser mejores. En el siguiente ejemplo, se muestra cómo automatizar la búsqueda.

Ajustar y comparar múltiples distribuciones

En el siguiente ejemplo se muestra cómo automatizar el ajuste y comparación de las múltiples distribuciones disponibles en scipy.stats. El código ha de permitir:

  • Ajustar todas las distribuciones disponibles en scipy.stats.

  • Poder preseleccionar un subconjunto de distribuciones candidatas en función de su soporte.

  • Mostrar los parámetros de cada ajuste.

  • Calcular los valores AIC y BIC para poder seleccionar la distribución con mejor ajuste.

  • Representación gráfica de los resultados

Funciones

Funciones empleadas para comparar múltiples distribuciones.

def seleccionar_distribuciones(
    familia: str = 'realall', 
    exclusiones: list = None,
    verbose: bool = True
) -> list:
    '''
    Esta función selecciona un subconjunto de las distribuciones disponibles
    en scipy.stats
    
    Parameters
    ----------
    familia : {'realall', 'realline', 'realplus', 'real0to1', 'discreta'}
        realall: distribuciones de la familia `realline` + `realplus`
        realline: distribuciones continuas en el soporte (-inf, +inf)
        realplus: distribuciones continuas en el soporte [0, +inf)
        real0to1: distribuciones continuas en el soporte [0,1]
        discreta: distribuciones discretas
        
    exclusiones : list, optional
        Lista de nombres de distribuciones a excluir. Por defecto se excluyen
        ['levy_stable', 'vonmises', 'erlang', 'studentized_range']
        
    verbose : bool
        Si se muestra la información de las distribuciones seleccionadas
        (the default `True`).
        
    Returns
    -------
    distribuciones: list
        listado con las distribuciones (los objetos) seleccionados.

    '''
    
    familias_validas = ['realall', 'realline', 'realplus', 'real0to1', 'discreta']
    if familia not in familias_validas:
        raise ValueError(f"El parámetro 'familia' debe ser uno de {familias_validas}, "
                         f"pero se recibió '{familia}'")
    
    # Obtener todas las distribuciones
    distribuciones = [getattr(stats, d) for d in dir(stats) 
                     if isinstance(getattr(stats, d), (stats.rv_continuous, stats.rv_discrete))]
    
    # Exclusiones por defecto si no se especifican
    if exclusiones is None:
        exclusiones = ['levy_stable', 'vonmises', 'erlang', 'studentized_range']
    
    exclusiones_set = set(exclusiones)
    distribuciones = [dist for dist in distribuciones if dist.name not in exclusiones_set]

    info_distribuciones = pd.DataFrame([
        {
            'distribucion': dist.name,
            'tipo': 'continua' if isinstance(dist, stats.rv_continuous) else 'discreta',
            'soporte_inf': dist.a,
            'soporte_sup': dist.b,
            'dist_obj': dist
        }
        for dist in distribuciones
    ]).sort_values(by=['soporte_inf', 'soporte_sup']).reset_index(drop=True)
    
    soportes = {
        'realline': [-np.inf, np.inf],
        'realplus': [0, np.inf],
        'real0to1': [0, 1]
    }
    
    if familia == 'realall':
        mask = (info_distribuciones['tipo'] == 'continua') & (
            ((info_distribuciones['soporte_inf'] == -np.inf) & 
             (info_distribuciones['soporte_sup'] == np.inf)) |
            ((info_distribuciones['soporte_inf'] == 0) & 
             (info_distribuciones['soporte_sup'] == np.inf))
        )
        info_distribuciones = info_distribuciones[mask].reset_index(drop=True)
        
    elif familia in soportes:
        mask = (info_distribuciones['tipo'] == 'continua') & \
               (info_distribuciones['soporte_inf'] == soportes[familia][0]) & \
               (info_distribuciones['soporte_sup'] == soportes[familia][1])
        info_distribuciones = info_distribuciones[mask].reset_index(drop=True)
        
    elif familia == 'discreta':
        info_distribuciones = info_distribuciones[
            info_distribuciones['tipo'] == 'discreta'
        ].reset_index(drop=True)
    
    seleccion = info_distribuciones['dist_obj'].tolist()
    info_distribuciones_display = info_distribuciones.drop(columns=['dist_obj'])
    
    if not seleccion:
        print("⚠️  No se encontraron distribuciones para la familia especificada")
    
    if verbose:
        print("---------------------------------------------------")
        print("       Distribuciones seleccionadas                ")
        print("---------------------------------------------------")
        with pd.option_context('display.max_rows', None, 'display.max_columns', None): 
            print(info_distribuciones_display)
    
    return seleccion


def comparar_distribuciones(
    x: np.ndarray,
    familia: str = 'realall', 
    criterio: str = 'aic',
    exclusiones: list = None,
    verbose: bool = True,
    suppress_warnings: bool = True
) -> pd.DataFrame:
    '''
    Esta función selecciona y ajusta un subconjunto de las distribuciones 
    disponibles en scipy.stats. Para cada distribución calcula los valores de
    Log Likelihood, AIC y BIC.
    
    Parameters
    ----------
    x : array_like
        datos con los que ajustar la distribución.
        
    familia : {'realall', 'realline', 'realplus', 'real0to1', 'discreta'}
        realall: distribuciones de la familia `realline` + `realplus`
        realline: distribuciones continuas en el soporte (-inf, +inf)
        realplus: distribuciones continuas en el soporte [0, +inf)
        real0to1: distribuciones continuas en el soporte [0,1]
        discreta: distribuciones discretas
    
    criterio : {'aic', 'bic'}
        criterio de ordenación de mejor a peor ajuste.
    
    exclusiones : list, optional
        Lista de nombres de distribuciones a excluir. Por defecto se excluyen
        ['levy_stable', 'vonmises', 'erlang', 'studentized_range']
    
    verbose : bool
        Si se muestra la información de las distribuciones seleccionadas
        (the default `True`).
    
    suppress_warnings : bool
        Si se suprimen los warnings numéricos durante el ajuste (RuntimeWarning).
        Se recomienda True para obtener salida limpia (the default `True`).
        
    Returns
    -------
    resultados: pd.DataFrame
        distribucion: nombre de la distribución.
        log_likelihood: logaritmo del likelihood del ajuste.
        aic: métrica AIC.
        bic: métrica BIC.
        n_parametros: número de parámetros de la distribución.
        parametros: parámetros tras el ajuste

    '''
    
    criterios_validos = ['aic', 'bic']
    if criterio not in criterios_validos:
        raise ValueError(f"El parámetro 'criterio' debe ser uno de {criterios_validos}, "
                         f"pero se recibió '{criterio}'")
    
    if exclusiones is not None and not isinstance(exclusiones, (list, tuple, set)):
        raise TypeError("El parámetro 'exclusiones' debe ser una lista, tupla o set")
    
    x = np.asarray(x)
    if x.size == 0:
        raise ValueError("El array 'x' no puede estar vacío")
    
    if not np.all(np.isfinite(x)):
        raise ValueError("El array 'x' contiene valores NaN o infinitos")
    
    distribuciones = seleccionar_distribuciones(
        familia=familia, 
        exclusiones=exclusiones, 
        verbose=False
    )
    
    resultados = []
    for distribucion in tqdm(distribuciones):
        
        if verbose:
            print(f"Ajustando distribución: {distribucion.name}")
        
        try:
            if suppress_warnings:
                with warnings.catch_warnings():
                    warnings.simplefilter("ignore")
                    parametros = distribucion.fit(data=x)
            else:
                parametros = distribucion.fit(data=x)
            
            # Determinar si es continua o discreta
            is_continuous = isinstance(distribucion, stats.rv_continuous)
            sig_func = distribucion._pdf if is_continuous else distribucion._pmf
            log_func = distribucion.logpdf if is_continuous else distribucion.logpmf
            
            # Extraer nombres de parámetros
            nombre_parametros = [p for p in inspect.signature(sig_func).parameters 
                               if p != 'x'] + ["loc", "scale"]
            parametros_dict = dict(zip(nombre_parametros, parametros))
            
            # Calcular métricas
            log_likelihood = log_func(x, *parametros).sum()
            aic = -2 * log_likelihood + 2 * len(parametros)
            bic = -2 * log_likelihood + np.log(x.shape[0]) * len(parametros)
            
            # Almacenar resultados
            resultados.append({
                'distribucion': distribucion.name,
                'log_likelihood': log_likelihood,
                'aic': aic,
                'bic': bic,
                'n_parametros': len(parametros),
                'parametros': parametros_dict
            })
            
        except (ValueError, RuntimeError, FloatingPointError) as e:
            if verbose:
                print(f"Error al ajustar {distribucion.name}: {type(e).__name__}: {e}\n")
        except Exception as e:
            if verbose:
                print(f"Error inesperado al ajustar {distribucion.name}: {type(e).__name__}: {e}\n")
    
    # Verificar si se ajustó alguna distribución
    if not resultados and verbose:
        print("\n⚠️  No se pudo ajustar ninguna distribución con éxito.")
        resultados = pd.DataFrame(
            columns=['distribucion', 'log_likelihood', 'aic', 'bic', 'n_parametros', 'parametros']
        )
    else:
        resultados = pd.DataFrame(resultados).sort_values(by=criterio).reset_index(drop=True)
    
    return resultados

Datos

De nuevo se emplean como datos el precio de los diamantes disponible en data set diamonds de la librería seaborn, en concreto, la columna price.

# Datos
# ==============================================================================
datos = sns.load_dataset('diamonds')
datos = datos.loc[datos.cut == 'Fair', 'price']

Ajuste de distribuciones

# Ajuste y comparación de distribuciones
# ==============================================================================
resultados = comparar_distribuciones(
                x=datos.to_numpy(),
                familia='realall',
                criterio='aic',
                verbose=False,
                suppress_warnings=True
            )
resultados

  0%|          | 0/84 [00:00<?, ?it/s]
distribucion log_likelihood aic bic n_parametros parametros
0 recipinvgauss -1.488310e+04 2.977220e+04 2.978835e+04 3 {'mu': 0.8564617364318987, 'loc': 75.589562124...
1 geninvgauss -1.488228e+04 2.977256e+04 2.979410e+04 4 {'p': 0.03895738587977514, 'b': 1.373372982791...
2 norminvgauss -1.488307e+04 2.977415e+04 2.979569e+04 4 {'a': 53.047412557163966, 'b': 53.025400276075...
3 genhyperbolic -1.488252e+04 2.977503e+04 2.980195e+04 5 {'p': 0.06748200212521904, 'a': 19.46831101851...
4 lognorm -1.488561e+04 2.977722e+04 2.979337e+04 3 {'s': 0.775749425283541, 'loc': 30.76091834294...
... ... ... ... ... ... ...
78 exponweib -1.697787e+04 3.396375e+04 3.398529e+04 4 {'a': 3.6435974274790723, 'c': 0.1273335852525...
79 invweibull -1.799302e+04 3.599205e+04 3.600820e+04 3 {'c': 0.1931244916343704, 'loc': 336.999783682...
80 ncx2 -1.024232e+06 2.048472e+06 2.048493e+06 4 {'df': 0.0011534145941760432, 'nc': 0.34024851...
81 truncnorm -inf inf inf 4 {'a': 388.84586460519677, 'b': 506.69039671197...
82 weibull_min -inf inf inf 3 {'c': 1.1567328971778637, 'loc': 495.039530173...

83 rows × 6 columns

Gráficos

def fit_plot_distribution(
    x: np.ndarray, 
    nombre_distribucion, 
    ax=None, 
    verbose: bool = True,
):
    '''
    Esta función superpone la(s) curva(s) de densidad de una o varias 
    distribuciones con el histograma de los datos.
    
    Parameters
    ----------
    x : array_like
        datos con los que ajustar la distribución.
        
    nombre_distribucion : str o list
        nombre de una distribución o lista de nombres de distribuciones 
        disponibles en `scipy.stats`.
    
    ax : matplotlib.axes.Axes, optional
        eje sobre el que realizar el gráfico. Si es None, se crea uno nuevo.
    
    verbose : bool, optional
        Si True, imprime los resultados del ajuste. Para múltiples distribuciones,
        imprime advertencias de errores (default True).
        
    Returns
    -------
    ax: matplotlib.axes.Axes
        gráfico creado
        
    '''
    
    # Validación de entrada
    x = np.asarray(x)
    if x.size == 0:
        raise ValueError("El array 'x' no puede estar vacío")
    if not np.all(np.isfinite(x)):
        raise ValueError("El array 'x' contiene valores NaN o infinitos")
    
    if isinstance(nombre_distribucion, str):
        nombre_distribucion = [nombre_distribucion]
    if not isinstance(nombre_distribucion, (list, tuple)):
        raise TypeError("El parámetro 'nombre_distribucion' debe ser un string, lista o tupla")
    
    # Crear gráfico
    if ax is None:
        fig, ax = plt.subplots(figsize=(7,4))
    
    # Graficar histograma
    ax.hist(x=x, density=True, bins=30, color="#3182bd", alpha=0.5)
    ax.plot(x, np.full_like(x, -0.01), '|k', markeredgewidth=1)
    ax.set_title('Ajuste distribución' if len(nombre_distribucion) == 1 else 
                 'Ajuste distribuciones')
    ax.set_xlabel('x')
    ax.set_ylabel('Densidad de probabilidad')
    
    for nombre in tqdm(nombre_distribucion):
        if not hasattr(stats, nombre):
            if verbose:
                print(f"⚠️  La distribución '{nombre}' no existe en scipy.stats")
            continue
        
        try:

            distribucion = getattr(stats, nombre)
            parametros = distribucion.fit(data=x)
            
            is_continuous = isinstance(distribucion, stats.rv_continuous)
            pdf_func = distribucion.pdf if is_continuous else distribucion.pmf
            
            x_hat = np.linspace(min(x), max(x), num=100)
            y_hat = pdf_func(x_hat, *parametros)
            ax.plot(x_hat, y_hat, linewidth=2, label=distribucion.name)
            
        except Exception as e:
            if verbose:
                print(f"⚠️  Error al ajustar la distribución '{nombre}': {type(e).__name__}: {e}")
    
        ax.legend(loc='best')

    return ax

Se muestra la mejor distribución acorde al criterio AIC.

fig, ax = plt.subplots(figsize=(6, 3))

fit_plot_distribution(
    x=datos.to_numpy(),
    nombre_distribucion=resultados['distribucion'][0],
    ax=ax
);

  0%|          | 0/1 [00:00<?, ?it/s]

Las curvas de densidad de probabilidad para las 5 mejores distribuciones.

fig, ax = plt.subplots(figsize=(6, 3))

fit_plot_distribution(
    x=datos.to_numpy(),
    nombre_distribucion=resultados['distribucion'][:5].tolist(),
    ax=ax
);

  0%|          | 0/5 [00:00<?, ?it/s]

Función de densidad, cuantil y muestreo

Todas las funciones implementadas en scipy.stats disponen de los métodos pdf(), logpdf(), cdf(), ppf() y rvs() con los que calcular la densidad, logaritmo de densidad, probabilidad acumulada, cuantiles, y muestreo de nuevos valores. Por ejemplo, se pueden simular 5 nuevos valores de diamantes acorde a la distribución johnsonsu.

# Definición de la distribución
distribucion = stats.johnsonsu

# Ajuste para obtener el valor de los parámetros
parametros   = distribucion.fit(datos.to_numpy())

# Muestreo aleatorio
distribucion.rvs(*parametros, size=5)

array([1716.35766217, 3479.75678645, 2604.21442503, 1218.40672023,
       2507.58702509])

Información de sesión

import session_info
session_info.show(html=False)

-----
matplotlib          3.10.8
numpy               2.2.6
pandas              2.3.3
scipy               1.15.3
seaborn             0.13.2
session_info        v1.0.1
statsmodels         0.14.6
tqdm                4.67.1
-----
IPython             9.8.0
jupyter_client      8.7.0
jupyter_core        5.9.1
-----
Python 3.13.11 | packaged by Anaconda, Inc. | (main, Dec 10 2025, 21:28:48) [GCC 14.3.0]
Linux-6.14.0-37-generic-x86_64-with-glibc2.39
-----
Session information updated at 2026-01-16 12:32

Bibliografía

Distributions for Modeling Location, Scale, and Shape Using GAMLSS in R By Robert A. Rigby, Mikis D. Stasinopoulos, Gillian Z. Heller, Fernanda De Bastiani.

Delignette-Muller, M., & Dutang, C. (2015). fitdistrplus: An R Package for Fitting Distributions. doi:http://dx.doi.org/10.18637/jss.v064.i04

Instrucciones para citar

¿Cómo citar este documento?

Si utilizas este documento o alguna parte de él, te agradecemos que lo cites. ¡Muchas gracias!

Bootstrapping para intervalos de confianza y contraste de hipótesis por Joaquín Amat Rodrigo, disponible bajo una licencia Attribution-NonCommercial-ShareAlike 4.0 International (CC BY-NC-SA 4.0 DEED) en https://www.cienciadedatos.net/documentos/pystats01-ajuste-distribuciones-python

¿Te ha gustado el artículo? Tu ayuda es importante

Tu contribución me ayudará a seguir generando contenido divulgativo gratuito. ¡Muchísimas gracias! 😊

Become a GitHub Sponsor

Creative Commons Licence

Este documento creado por Joaquín Amat Rodrigo tiene licencia Attribution-NonCommercial-ShareAlike 4.0 International.

Se permite:

  • Compartir: copiar y redistribuir el material en cualquier medio o formato.

  • Adaptar: remezclar, transformar y crear a partir del material.

Bajo los siguientes términos:

  • Atribución: Debes otorgar el crédito adecuado, proporcionar un enlace a la licencia e indicar si se realizaron cambios. Puedes hacerlo de cualquier manera razonable, pero no de una forma que sugiera que el licenciante te respalda o respalda tu uso.

  • No-Comercial: No puedes utilizar el material para fines comerciales.

  • Compartir-Igual: Si remezclas, transformas o creas a partir del material, debes distribuir tus contribuciones bajo la misma licencia que el original.