2010-07-13 13 views
7

Ich habe einen Datensatz, von dem ich weiß, dass er eine Pareto-Verteilung hat. Kann mir jemand zeigen, wie man diesen Datensatz in Scipy einpasst? Ich habe den folgenden Code zum laufen, aber ich habe keine Ahnung, was mir zurückgegeben wird (a, b, c). Nach dem Erhalt von a, b, c, wie berechne ich dann die Varianz mit ihnen?Anpassen einer Pareto-Verteilung mit (Python) Scipy

Antwort

1

Seien Sie sehr vorsichtig bei der Anpassung der Stromgesetze !! Viele gemeldete Machtgesetze sind durch ein Potenzgesetz tatsächlich schlecht angepasst. Alle Details finden Sie unter Clauset et al. (auch unter arxiv, wenn Sie keinen Zugriff auf das Journal haben). Sie haben einen companion website zu dem Artikel, der jetzt eine Verbindung zu einer Python-Implementierung herstellt. Ich weiß nicht, ob es Scipy benutzt, weil ich die R-Implementierung benutzt habe, als ich es zuletzt benutzt habe.

+1

Die Python-Implementierung (http://code.google.com/p/agpy/wiki/PowerLaw) enthält zwei Versionen; einer hängt von der Zahl ab, einer nicht. (Ich schrieb es) – keflavich

3

Hier ist eine schnell geschriebene Version, die einige Hinweise von der Referenzseite nimmt, die Rupert gab. Dies ist derzeit in scipy und statsmodels in Arbeit und erfordert MLE mit einigen festen oder eingefrorenen Parametern, die nur in den Stammversionen verfügbar sind. Es sind noch keine Standardfehler für die Parameterschätzer oder andere Ergebnisstatistiken verfügbar.

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