2010-07-13 10 views
7

Ho un set di dati che conosco ha una distribuzione Pareto. Qualcuno può indicarmi come montare questo set di dati in Scipy? Ho ottenuto il codice seguente per eseguire ma non ho idea di cosa mi viene restituito (a, b, c). Inoltre, dopo aver ottenuto a, b, c, come faccio a calcolare la varianza che li usa?Inserimento di una distribuzione di pareto con (python) Scipy

import scipy.stats as ss 
import scipy as sp 

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

risposta

1

Fare molta attenzione alle leggi di alimentazione !! Molte leggi di potere riportate sono in realtà malamente montate da una legge di potenza. Vedi Clauset et al. per tutti i dettagli (anche su arxiv se non si ha accesso al giornale). Hanno un numero companion website per l'articolo che ora si collega a un'implementazione Python. Non so se usa Scipy perché ho usato la loro implementazione R quando l'ho usato per l'ultima volta.

+1

L'implementazione di Python (http://code.google.com/p/agpy/wiki/PowerLaw) include due versioni; uno dipende da Numpy, uno no. (L'ho scritto io) – keflavich

3

Ecco una versione scritta rapidamente, prendendo alcuni suggerimenti dalla pagina di riferimento che Rupert ha dato. Questo è attualmente in lavorazione in scipy e statsmodels e richiede MLE con alcuni parametri fissi o congelati, che è disponibile solo nelle versioni trunk. Non sono ancora disponibili errori standard sulle stime dei parametri o altre statistiche sui risultati.

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