Ajuste de una distribución pareto con (python) Scipy

Question

Tengo un conjunto de datos que sé que tiene una distribución de Pareto. ¿Alguien puede indicarme cómo ajustar este conjunto de datos en Scipy? Tengo el siguiente código para ejecutar, pero no tengo idea de lo que me está devolviendo (a, b, c). Además, después de obtener a, b, c, ¿cómo calculo la varianza al usarlos?

import scipy.stats as ss 
import scipy as sp 

a,b,c=ss.pareto.fit(data) 

Answer 1

Sea muy cuidadoso al ajustar las leyes de potencia !! Muchas leyes de energía informadas en realidad están mal ajustadas por una ley de energía. Consulte Clauset et al. para obtener todos los detalles (también en arxiv si no tiene acceso al diario). Tienen un companion website para el artículo que ahora se vincula a una implementación de Python. No sé si usa Scipy porque utilicé su implementación R la última vez que lo usé.

Answer 2

Aquí hay una versión escrita rápidamente, que toma algunas sugerencias de la página de referencia que dio Rupert. Esto es actualmente un trabajo en progreso en scipy y modelos de estadísticas y requiere MLE con algunos parámetros fijos o congelados, que solo está disponible en las versiones troncales. Aún no hay errores estándar en las estimaciones de parámetros u otras estadísticas de resultados.

'''estimating pareto with 3 parameters (shape, loc, scale) with nested 
minimization, MLE inside minimizing Kolmogorov-Smirnov statistic 

running some examples looks good 
Author: josef-pktd 
''' 

import numpy as np 
from scipy import stats, optimize 
#the following adds my frozen fit method to the distributions 
#scipy trunk also has a fit method with some parameters fixed. 
import scikits.statsmodels.sandbox.stats.distributions_patch 

true = (0.5, 10, 1.) # try different values 
shape, loc, scale = true 
rvs = stats.pareto.rvs(shape, loc=loc, scale=scale, size=1000) 

rvsmin = rvs.min() #for starting value to fmin 


def pareto_ks(loc, rvs): 
    est = stats.pareto.fit_fr(rvs, 1., frozen=[np.nan, loc, np.nan]) 
    args = (est[0], loc, est[1]) 
    return stats.kstest(rvs,'pareto',args)[0] 

locest = optimize.fmin(pareto_ks, rvsmin*0.7, (rvs,)) 
est = stats.pareto.fit_fr(rvs, 1., frozen=[np.nan, locest, np.nan]) 
args = (est[0], locest[0], est[1]) 
print 'estimate' 
print args 
print 'kstest' 
print stats.kstest(rvs,'pareto',args) 
print 'estimation error', args - np.array(true)