2010-07-13 115 views
7

我有一个数据集,我知道帕累托分布。有人能指出我如何在Scipy中适合这个数据集吗?我得到了下面的代码运行,但我不知道什么是返回给我(a,b,c)。另外,获得a,b,c后,如何使用它们计算方差?拟合帕累托分布(python)Scipy

import scipy.stats as ss 
import scipy as sp 

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

回答

1

要非常小心的适合电源法则!许多报道的权力法律实际上严格遵守幂律。所有细节请参阅Clauset et al.(如果您无权访问期刊,也请登录arxiv)。对于现在链接到Python实现的文章,他们有companion website。不知道它是否使用Scipy,因为我上次使用它时使用了它的R实现。

+1

python实现(http://code.google.com/p/agpy/wiki/PowerLaw)包含两个版本;一个取决于numpy,一个不取决于。 (我写的) – keflavich 2012-01-12 00:05:58

3

下面是一个快速写入的版本,从鲁珀特给出的参考页面提供了一些提示。 这是目前正在进行中的scipy和statsmodels,并要求MLE的一些固定或冻结参数,它只在中继版本中可用。 参数估算值或其他结果统计数据尚未提供标准误差。

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