Ajuste y selección de distribuciones con Python

Joaquín Amat Rodrigo
Enero, 2020

Introducción


Identificar el tipo de distribución que tiene a una variable es un paso fundamental en prácticamente todos los estudios que implican datos, desde los contrastes de hipótesis hasta la creación de modelos por aprendizaje estadístico y machine learning.

En Python existen varias librerías que permiten ajustar distribuciones. En este documento se muestran las funcionalidades del módulo scipy.stats, haciendo hincapié en cómo comparar múltiples distribuciones con el objetivo de identificar a cuál de ellas se ajustan mejor los datos.



Ajuste de una distribución


Ajustar una distribución paramétrica a partir de un conjunto de datos consiste en encontrar el valor de los parámetros con los que, con mayor probabilidad, dicha distribución puede haber generado los datos observados. Por ejemplo, la distribución normal tiene dos parámetros (media y varianza), una vez conocidos estos dos parámetros, se conoce toda la distribución.

Existen varios métodos que permiten encontrar los parámetros óptimos que mejor se ajustan a los datos, uno de los más utilizados, y el que implementa scipy.stats, es el método de Maximum Likelihood Estimation (MLE) (máxima verosimilitud). scipy.stats dispone de más de 90 distribuciones, puede encontrarse un listado de todas ellas en:

In [50]:
from scipy import stats
import pandas as pd

Además de diferenciar entre distribuciones continuas y discretas, es útil poder seleccionarlas por el rango de valores sobre el que está definida cada distribución (dominio). Por ejemplo, si se quiere modelar la velocidad del viento, aunque no se conozca el tipo exacto de distribución, se puede acotar a aquellas cuyo rango de valores está limitado entre $0$ y $+\inf$.

In [51]:
# Distribuciones agrupadas por dominio
# ==============================================================================
distribuciones = [getattr(stats,d) for d in dir(stats) \
                  if isinstance(getattr(stats,d), (stats.rv_continuous, stats.rv_discrete))]

distribucion = []     
dominio_1 = []
dominio_2 = []

for dist in distribuciones:
    distribucion.append(dist.name)
    dominio_1.append(dist.a)
    dominio_2.append(dist.b)
    
info_distribuciones = pd.DataFrame({
                        'distribucion': distribucion,
                        'dominio_1': dominio_1,
                        'dominio_2': dominio_2
                      })

info_distribuciones = info_distribuciones \
                      .sort_values(by=['dominio_1', 'dominio_2'])\
                      .reset_index(drop=True)

print("-------------------------------------")
print("Información distribuciones scipy.stat")
print("-------------------------------------")
display(info_distribuciones)
-------------------------------------
Información distribuciones scipy.stat
-------------------------------------
distribucion dominio_1 dominio_2
0 frechet_l -inf 0.0
1 levy_l -inf 0.0
2 weibull_max -inf 0.0
3 cauchy -inf inf
4 crystalball -inf inf
... ... ... ...
111 geom 1.0 inf
112 logser 1.0 inf
113 pareto 1.0 inf
114 yulesimon 1.0 inf
115 zipf 1.0 inf

116 rows × 3 columns

Comparación de ajustes


El método de ajuste por Maximum Likelihood Estimation (MLE) consigue, dada una determinada distribución, encontrar el valor sus parámetros con los que mejor se ajusta a los datos. En la práctica, es frecuente no conocer de antemano qué tipo de distribución siguen los datos. Aunque con la experiencia, el analista suele poder acotar las distribuciones candidatas (distribuciones continuas o discretas, distribuciones solo positivas, etc), es necesaria una forma de cuantificar la bondad de ajuste de cada distribución y, además, poder comparar entre distintas distribuciones candidatas.

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$.

Para todas ellas, cuanto menor sea el valor, mejor el ajuste. El cambio en la dirección de selección respecto al log likelihood se debe a que en ambas métricas se multiplica por -2, por lo tanto, si interesan valores elevados de log likelihood, al multiplicar por -2, interesan valores muy negativos.

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 1


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:

In [52]:
# Tratamiento de datos
# ==============================================================================
import pandas as pd
import numpy as np
import seaborn as sns

# Gráficos
# ==============================================================================
import matplotlib.pyplot as plt
from matplotlib import style

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

# Configuración matplotlib
# ==============================================================================
#plt.rcParams['image.cmap'] = "bwr"
#plt.rcParams['figure.dpi'] = "100"
plt.rcParams['savefig.bbox'] = "tight"
style.use('ggplot') or plt.style.use('ggplot')

# Configuración warnings
# ==============================================================================
import warnings
warnings.filterwarnings('ignore')

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.

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

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().

In [54]:
# Estadísticos descriptivos
# ==============================================================================
datos.describe()
Out[54]:
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
In [55]:
# Gráficos distribución observada (empírica)
# ==============================================================================
fig, axs = plt.subplots(nrows=1, ncols=2, figsize=(10, 4))

# 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();
In [56]:
# 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=(7,4))
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"Dominio:        {[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
Dominio:        [-inf, inf]
Parámetros:     {'loc': 4358.757763975155, 'scale': 3559.2807303891086}
Log likelihood: -15449.966194325283
AIC:            30903.932388650566
BIC:            30914.700367566522

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

In [57]:
# Ajuste distribución normal
#===============================================================================
# 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=(7,4))
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"Dominio:        {[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
Dominio:        [0.0, inf]
Parámetros:     {'a': 14.399103124683087, 'loc': -6244.786483560789, 'scale': 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 es 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 dominio.

  • 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

Código

Funciones empleadas para comparar múltiples distribuciones.

In [58]:
from scipy import stats
import pandas as pd
import numpy as np
import tqdm
import inspect
import warnings
warnings.filterwarnings('ignore')

def seleccionar_distribuciones(familia='realall', verbose=True):
    '''
    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 dominio (-inf, +inf)
        realplus: distribuciones continuas en el dominio [0, +inf)
        real0to1: distribuciones continuas en el dominio [0,1]
        discreta: distribuciones discretas
        
    verbose : bool
        Si se muestra información de las distribuciones seleccionadas
        (the default `True`).
        
    Returns
    -------
    distribuciones: list
        listado con las distribuciones (los objetos) seleccionados.
        
    Raises
    ------
    Exception
        Si `familia` es distinto de 'realall', 'realline', 'realplus', 'real0to1',
        o 'discreta'.
        
    Notes
    -----
        Las distribuciones levy_stable y vonmises han sido excluidas por el momento.

    '''
    
    distribuciones = [getattr(stats,d) for d in dir(stats) \
                     if isinstance(getattr(stats,d), (stats.rv_continuous, stats.rv_discrete))]
    
    exclusiones = ['levy_stable', 'vonmises']
    distribuciones = [dist for dist in distribuciones if dist.name not in exclusiones]
            
    dominios = {
        'realall' : [-np.inf, np.inf],
        'realline': [np.inf,np.inf],
        'realplus': [0, np.inf],
        'real0to1': [0, 1], 
        'discreta': [None, None],
    }

    distribucion = []
    tipo = []
    dominio_inf = []
    dominio_sup = []

    for dist in distribuciones:
        distribucion.append(dist.name)
        tipo.append(np.where(isinstance(dist, stats.rv_continuous), 'continua', 'discreta'))
        dominio_inf.append(dist.a)
        dominio_sup.append(dist.b)
    
    info_distribuciones = pd.DataFrame({
                            'distribucion': distribucion,
                            'tipo': tipo,
                            'dominio_inf': dominio_inf,
                            'dominio_sup': dominio_sup
                          })

    info_distribuciones = info_distribuciones \
                          .sort_values(by=['dominio_inf', 'dominio_sup'])\
                          .reset_index(drop=True)
    
    if familia in ['realall', 'realline', 'realplus', 'real0to1']:
        info_distribuciones = info_distribuciones[info_distribuciones['tipo']=='continua']
        condicion = (info_distribuciones['dominio_inf'] == dominios[familia][0]) & \
                    (info_distribuciones['dominio_sup'] == dominios[familia][1]) 
        info_distribuciones = info_distribuciones[condicion].reset_index(drop=True)
        
    if familia in ['discreta']:
        info_distribuciones = info_distribuciones[info_distribuciones['tipo']=='discreta']
        
    seleccion = [dist for dist in distribuciones \
                 if dist.name in info_distribuciones['distribucion'].values]
    
    
    if verbose:
        print("---------------------------------------------------")
        print("       Distribuciones seleccionadas                ")
        print("---------------------------------------------------")
        with pd.option_context('display.max_rows', None, 'display.max_columns', None): 
            print(info_distribuciones)
    
    return seleccion


def comparar_distribuciones(x, familia='realall', ordenar='aic', verbose=True):
    '''
    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 dominio (-inf, +inf)
        realplus: distribuciones continuas en el dominio [0, +inf)
        real0to1: distribuciones continuas en el dominio [0,1]
        discreta: distribuciones discretas
    
    ordenar : {'aic', 'bic'}
        criterio de ordenación de mejor a peor ajuste.
    
    verbose : bool
        Si se muestra información de las distribuciones seleccionadas
        (the default `True`).
        
    Returns
    -------
    resultados: data.frame
        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 de la distribución.
        parametros: parámetros del tras el ajuste
        
    Raises
    ------
    Exception
        Si `familia` es distinto de 'realall', 'realline', 'realplus', 'real0to1',
        o 'discreta'.
        
    Notes
    -----

    '''
    
    distribuciones = seleccionar_distribuciones(familia=familia, verbose=verbose)
    distribucion_ = []
    log_likelihood_= []
    aic_ = []
    bic_ = []
    n_parametros_ = []
    parametros_ = []
    
    for i, distribucion in enumerate(distribuciones):
        
        print(f"{i+1}/{len(distribuciones)} Ajustando distribución: {distribucion.name}")
        
        try:
            parametros = distribucion.fit(data=x)
            nombre_parametros = [p for p in inspect.signature(distribucion._pdf).parameters \
                                 if not p=='x'] + ["loc","scale"]
            parametros_dict = dict(zip(nombre_parametros, parametros))
            log_likelihood = distribucion.logpdf(x, *parametros).sum()
            aic = -2 * log_likelihood + 2 * len(parametros)
            bic = -2 * log_likelihood + np.log(x.shape[0]) * len(parametros)
            
            distribucion_.append(distribucion.name)
            log_likelihood_.append(log_likelihood)
            aic_.append(aic)
            bic_.append(bic)
            n_parametros_.append(len(parametros))
            parametros_.append(parametros_dict)
            
            resultados = pd.DataFrame({
                            'distribucion': distribucion_,
                            'log_likelihood': log_likelihood_,
                            'aic': aic_,
                            'bic': bic_,
                            'n_parametros': n_parametros_,
                            'parametros': parametros_,
                
                         })
            
            resultados = resultados.sort_values(by=ordenar).reset_index(drop=True)
            
        except Exception as e:
            print(f"Error al tratar de ajustar la distribución {distribucion.name}")
            print(e)
            print("")
            
    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.

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

Ajuste de distribuciones

In [60]:
# Ajuste y comparación de distribuciones
# ==============================================================================
resultados = comparar_distribuciones(
                x=datos.to_numpy(),
                familia='realall',
                ordenar='aic',
                verbose=False
            )
resultados
1/28 Ajustando distribución: cauchy
2/28 Ajustando distribución: crystalball
3/28 Ajustando distribución: dgamma
4/28 Ajustando distribución: dweibull
5/28 Ajustando distribución: exponnorm
6/28 Ajustando distribución: genextreme
7/28 Ajustando distribución: genlogistic
8/28 Ajustando distribución: gennorm
9/28 Ajustando distribución: gumbel_l
10/28 Ajustando distribución: gumbel_r
11/28 Ajustando distribución: hypsecant
12/28 Ajustando distribución: johnsonsu
13/28 Ajustando distribución: kappa4
14/28 Ajustando distribución: laplace
15/28 Ajustando distribución: loggamma
16/28 Ajustando distribución: logistic
17/28 Ajustando distribución: loguniform
18/28 Ajustando distribución: moyal
19/28 Ajustando distribución: nct
20/28 Ajustando distribución: norm
21/28 Ajustando distribución: norminvgauss
22/28 Ajustando distribución: pearson3
23/28 Ajustando distribución: powernorm
24/28 Ajustando distribución: reciprocal
25/28 Ajustando distribución: skewnorm
26/28 Ajustando distribución: t
27/28 Ajustando distribución: truncnorm
28/28 Ajustando distribución: tukeylambda
Out[60]:
distribucion log_likelihood aic bic n_parametros parametros
0 johnsonsu -1.488566e+04 2.977932e+04 2.980085e+04 4 {'a': -6.072232355774776, 'b': 1.2885984469584...
1 norminvgauss -1.488888e+04 2.978576e+04 2.980729e+04 4 {'a': 167.4132134376864, 'b': 167.408351564553...
2 exponnorm -1.489309e+04 2.979219e+04 2.980834e+04 3 {'K': 15.41214099000929, 'loc': 743.8580247955...
3 pearson3 -1.490157e+04 2.980914e+04 2.982529e+04 3 {'skew': 1.6311412759431207, 'loc': 4358.75780...
4 genextreme -1.490368e+04 2.981335e+04 2.982950e+04 3 {'c': -0.3524573977513747, 'loc': 2547.4073584...
5 nct -1.493197e+04 2.987195e+04 2.989348e+04 4 {'df': 2.090423595048227, 'nc': 2.592582858699...
6 moyal -1.496393e+04 2.993187e+04 2.994264e+04 2 {'loc': 2529.6661504971357, 'scale': 1307.9215...
7 skewnorm -1.497556e+04 2.995711e+04 2.997326e+04 3 {'a': 38.123603745838736, 'loc': 553.312574042...
8 gumbel_r -1.504785e+04 3.009970e+04 3.011047e+04 2 {'loc': 2920.8689835486775, 'scale': 2178.9023...
9 genlogistic -1.504806e+04 3.010211e+04 3.011826e+04 3 {'c': 1462.8065290816894, 'loc': -12967.822516...
10 kappa4 -1.508691e+04 3.018182e+04 3.020335e+04 4 {'h': 1.474013806660743, 'k': 0.36166011975020...
11 powernorm -1.518749e+04 3.038099e+04 3.039714e+04 3 {'c': 1.4325348218947176e-05, 'loc': -989.9534...
12 tukeylambda -1.519527e+04 3.039654e+04 3.041270e+04 3 {'lam': -0.39734231519118934, 'loc': 3194.6385...
13 t -1.519599e+04 3.039799e+04 3.041414e+04 3 {'df': 1.9274620981470725, 'loc': 3207.5252784...
14 gennorm -1.522856e+04 3.046313e+04 3.047928e+04 3 {'beta': 0.778851883696015, 'loc': 3145.000000...
15 dweibull -1.523487e+04 3.047574e+04 3.049189e+04 3 {'c': 0.9252459732991334, 'loc': 3172.00000000...
16 dgamma -1.523747e+04 3.048095e+04 3.049710e+04 3 {'a': 0.9231119684247481, 'loc': 3144.99999999...
17 laplace -1.524740e+04 3.049879e+04 3.050956e+04 2 {'loc': 3282.0, 'scale': 2385.7888198757764}
18 cauchy -1.526197e+04 3.052795e+04 3.053871e+04 2 {'loc': 3024.3671767425167, 'scale': 1365.1725...
19 hypsecant -1.527175e+04 3.054750e+04 3.055827e+04 2 {'loc': 3557.045200428866, 'scale': 2030.08240...
20 logistic -1.531691e+04 3.063782e+04 3.064859e+04 2 {'loc': 3757.3286130383717, 'scale': 1750.5176...
21 norm -1.544997e+04 3.090393e+04 3.091470e+04 2 {'loc': 4358.757763975155, 'scale': 3559.28073...
22 loggamma -1.546426e+04 3.093453e+04 3.095068e+04 3 {'c': 1127.6629399668445, 'loc': -846013.89540...
23 gumbel_l -1.594097e+04 3.188594e+04 3.189670e+04 2 {'loc': 6392.790264700308, 'scale': 4809.11678...
24 crystalball -1.619239e+04 3.239278e+04 3.241432e+04 4 {'beta': 5.214204453514922, 'm': 2.00475959603...
25 truncnorm -inf inf inf 4 {'a': 0.9015482215897133, 'b': 1.6009849280375...
26 reciprocal -inf inf inf 4 {'a': 0.9015482215897133, 'b': 1.6009849280375...
27 loguniform -inf inf inf 4 {'a': 0.9015482215897133, 'b': 1.6009849280375...

Gráficos

In [61]:
def plot_distribucion(x, nombre_distribucion, ax=None):
    '''
    Esta función superpone la curva de densidad de una distribución con el
    histograma de los datos.
    
    Parameters
    ----------
    x : array_like
        datos con los que ajustar la distribución.
        
    nombre_distribuciones : str
        nombre de una de las distribuciones disponibles en `scipy.stats`.
        
    Returns
    -------
    resultados: matplotlib.ax
        gráfico creado
        
    Raises
    ------
        
    Notes
    -----
    '''

    distribucion = getattr(stats, nombre_distribucion)

    parametros = distribucion.fit(data=x)

    nombre_parametros = [p for p in inspect.signature(distribucion._pdf).parameters \
                         if not p=='x'] + ["loc","scale"]
    parametros_dict = dict(zip(nombre_parametros, parametros))

    log_likelihood = distribucion.logpdf(x, *parametros).sum()

    aic = -2 * log_likelihood + 2 * len(parametros)
    bic = -2 * log_likelihood + np.log(x.shape[0]) * len(parametros)

    x_hat = np.linspace(min(x), max(x), num=100)
    y_hat = distribucion.pdf(x_hat, *parametros)
    
    if ax is None:
        fig, ax = plt.subplots(figsize=(7,4))
        
    ax.plot(x_hat, y_hat, linewidth=2, label=distribucion.name)
    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')
    ax.set_xlabel('x')
    ax.set_ylabel('Densidad de probabilidad')
    ax.legend();

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


def plot_multiple_distribuciones(x, nombre_distribuciones, ax=None):
    '''
    Esta función superpone las curvas de densidad de varias distribuciones
    con el histograma de los datos.
    
    Parameters
    ----------
    x : array_like
        datos con los que ajustar la distribución.
        
    nombre_distribuciones : list
        lista con nombres de distribuciones disponibles en `scipy.stats`.
        
    Returns
    -------
    resultados: matplotlib.ax
        gráfico creado
        
    Raises
    ------
        
    Notes
    -----
    '''
    
    if ax is None:
        fig, ax = plt.subplots(figsize=(7,4))
        
    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 distribuciones')
    ax.set_xlabel('x')
    ax.set_ylabel('Densidad de probabilidad')
    
    for nombre in nombre_distribuciones:
        
        distribucion = getattr(stats, nombre)

        parametros = distribucion.fit(data=x)

        nombre_parametros = [p for p in inspect.signature(distribucion._pdf).parameters \
                             if not p=='x'] + ["loc","scale"]
        parametros_dict = dict(zip(nombre_parametros, parametros))

        log_likelihood = distribucion.logpdf(x, *parametros).sum()

        aic = -2 * log_likelihood + 2 * len(parametros)
        bic = -2 * log_likelihood + np.log(x.shape[0]) * len(parametros)

        x_hat = np.linspace(min(x), max(x), num=100)
        y_hat = distribucion.pdf(x_hat, *parametros)
        ax.plot(x_hat, y_hat, linewidth=2, label=distribucion.name)
    
    ax.legend();
    
    return ax

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

In [62]:
fig, ax = plt.subplots(figsize=(8,5))

plot_distribucion(
    x=datos.to_numpy(),
    nombre_distribucion=resultados['distribucion'][0],
    ax=ax
);
---------------------
Resultados del ajuste
---------------------
Distribución:   johnsonsu
Dominio:        [-inf, inf]
Parámetros:     {'a': -6.072232355774776, 'b': 1.2885984469584648, 'loc': 33.24040645468252, 'scale': 57.98885753098351}
Log likelihood: -14885.658863593013
AIC:            29779.317727186026
BIC:            29800.85368501794

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

In [63]:
fig, ax = plt.subplots(figsize=(8,5))

plot_multiple_distribuciones(
    x=datos.to_numpy(),
    nombre_distribuciones=resultados['distribucion'][:5],
    ax=ax
);

Resultados


Acorde al criterio AIC, las dos distribuciones que mejor se adaptan a los datos son: johnsonsu y norminvgauss.



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.

In [64]:
# 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)
Out[64]:
array([2791.58975437, 2100.10189361, 1866.55384315, 5681.39126389,
       2718.631243  ])

Información de sesión

In [65]:
from sinfo import sinfo
sinfo()
-----
ipykernel   5.3.4
matplotlib  3.3.2
numpy       1.19.2
pandas      1.1.3
scipy       1.5.2
seaborn     0.11.0
sinfo       0.3.1
statsmodels 0.12.1
tqdm        4.54.1
-----
IPython             7.18.1
jupyter_client      6.1.7
jupyter_core        4.6.3
jupyterlab          2.2.9
notebook            6.1.4
-----
Python 3.7.9 (default, Aug 31 2020, 12:42:55) [GCC 7.3.0]
Linux-5.4.0-1034-aws-x86_64-with-debian-buster-sid
4 logical CPU cores, x86_64
-----
Session information updated at 2021-01-07 12:01

Bibliografía


fitter

Statistics (scipy.stats)

univariateML

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

¿Cómo citar este documento?

Ajuste y selección de distribuciones con Python by Joaquín Amat Rodrigo, available under a Attribution 4.0 International (CC BY 4.0) at https://www.cienciadedatos.net/documentos/pystats01-ajuste-distribuciones-python.html


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

Mantener un sitio web tiene unos costes elevados, tu contribución me ayudará a seguir generando contenido divulgativo gratuito. ¡Muchísimas gracias! 😊


Creative Commons Licence
This work by Joaquín Amat Rodrigo is licensed under a Creative Commons Attribution 4.0 International License.