2010-07-13 20 views
7

J'ai un jeu de données dont je sais qu'il a une distribution Pareto. Quelqu'un peut-il me montrer comment adapter cet ensemble de données dans Scipy? J'ai le code ci-dessous pour fonctionner mais je n'ai aucune idée de ce qui me revient (a, b, c). De plus, après avoir obtenu a, b, c, comment puis-je calculer la variance en les utilisant?Ajustement d'une distribution pareto avec (python) Scipy

import scipy.stats as ss 
import scipy as sp 

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

Répondre

1

Soyez très prudent en ajustant les lois de puissance !! De nombreuses lois de puissance rapportées sont mal adaptées par une loi de puissance. Voir Clauset et al. pour tous les détails (également sur arxiv si vous n'avez pas accès au journal). Ils ont un companion website à l'article qui lie maintenant à une implémentation Python. Je ne sais pas si elle utilise Scipy parce que j'ai utilisé leur implémentation R quand je l'ai utilisé pour la dernière fois.

+1

L'implémentation Python (http://code.google.com/p/agpy/wiki/PowerLaw) inclut deux versions; on dépend de numpy, on ne le fait pas. (Je l'ai écrit) – keflavich

3

Voici une version rapidement écrite, en prenant quelques indications de la page de référence que Rupert a donné. Ceci est actuellement en cours dans scipy et statsmodels et nécessite MLE avec certains paramètres fixes ou gelés, qui est uniquement disponible dans les versions de tronc. Aucune erreur-type sur les estimations de paramètres ou autres statistiques de résultat n'est disponible pour le moment.

'''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)