參考文獻
Clauset, A., C.R. Shalizi, and M.E.J. Newman, Power-Law Distributions in Empirical Data. Siam Review, 2009. 51(4): p. 661-703.
Alstott, J., E.T. Bullmore, and D. Plenz, powerlaw: A Python Package for Analysis of Heavy-Tailed Distributions. Plos One, 2014. 9(1).
根據文獻一,對冪律分布的擬合主要包括:
1. 假定原始數據服從冪律分布,常見的方法通過直方圖,利用公式p(x) ∝ x-α,可以得知冪律分布服ln p(x) = αln x + constant, 其在雙對數坐標的圖像上會表現為一條直線,以此來初步觀測數據是否可能服從冪律分布。
2. 接著需要對Xmin和α進行估計
(α)【對於離散數據】
利用極大近似(maximum likelihood)的方法,能夠在有限的大樣本中得到比較精確的參數估計值,假定原始數據是從服從冪律分布的一組數據中截取X>Xmin的部分,我們可以利用極大似然估計(maximum likelihood estimators)對連續或者非連續的數據進行參數估計(α)【對於離散數據】
需要注意的是該方法並不會提示我們擬合結果是錯誤的,僅僅會返回對於數據的最佳擬合結果,並不會判斷冪律分布模型是否是對數據最佳的擬合模型,我們還需要後續的擬合檢驗。(同時建議數據應當保證n>50,可以得出精度較高的α)
對於Xmin
針對Xmin的估計採用的是Clauset et al的方法,可以同時應用於離散型和連續型數據,其方法的基本思想是選擇高於使得被測量數據的概率分布與最佳擬合冪律模型儘可能近似的Xmin估計。Xmin過大會大幅度的縮減原始數據集,有可能造成數據概率分布與冪律分布之間吻合度較差,另外Xmin過小會因為源數據與模型之間產生根本性的差異,使得分布發生變化,而在這中便是我們需要的Xmin估計。
3. 對數據下限估計結果的檢驗
4. 對冪律分布假設的檢驗
對冪律分布的假設檢驗主要是通過構建許多真正來自冪律分布的樣本數據,測量這些樣本數據與觀察數據之間的距離,如果兩者距離相差過大,那麼觀察數據就不能很好的冪律分布模型所擬合(not a plausible fit)。這裡需要注意採用何種方法來測量距離以及這種檢驗方法還是存在錯誤的可能,但是可以通過增加觀察數據容量,提高檢驗的精度。在距離的測量上,我們採用的是KS統計量,能夠產生較好的結果。但是對距離的顯著性,我們還需要進一步的計算。
標準的做法是用冪律模型對觀測數據進行擬合,並計算此時的KS統計量,接著基於最優擬合的參數(Xmin,α)構建冪律分布,以此產生足夠大的人工樣本數量,進而單獨的對人工樣本分布與觀測數據的冪律分布模型,計算它們的KS統計量,然後,我們簡單地計算出最終統計量大於經驗數據值的時間比例,這個比例就是p-value(Then we simply count what fraction of the time the resultingstatistic is larger than the value for the empirical data. This fraction is ourp-value)。需要注意的是,我們這裡計算的KS統計量,是人工樣本相較於原始數據的最佳冪律擬合結果,而不是數據的原始分布,也就是說用於KS統計量計算的是從原始數據中符合最佳擬合的數據子集。(we compute the KS statisticrelative to the best-fit power law for that data set, not relative to theoriginal distribution from which the data set was drawn.)
在我們的計算中,我們做出了一個相對保守的選擇:如果p <= 0,則冪律被排除。也就是說,如果百分之十或更少的概率,使我們僅偶然地獲得與模型一樣差的數據,就將其排除。在其他情況下,許多作者使用更優的規則p <= 0.05,但我們認為這會使某些候選分布真正通過冪律的可能性很小。當然,實際上,所採用的特定規則必須取決於研究者的判斷和當前的情況。In our calculations we have made the relatively conservative choicethat the power law is ruled out if p<=0.1that is, it is ruled out ifthere is a probability of 1 in 10 or less that we would merely by chance getdata that agree as poorly with the model as the data we have. (In othercontexts, many authors use the more lenient rule p<=0.05, but we feel this would letthrough some candidate distributions that have only a very small chance of reallyfollowing a power law. Of course, in practice, the particular rule adopted mustdepend on the judgment of the investigator and the circumstances at hand)
5. Direct comparison of models
如果通過了假設檢驗,我們還需要考慮,是否存在其他更好的分布模型能生成更好的擬合結果。
在這種情況下,存在可以直接比較兩個分布的方法,並且比KS測試更容易實現。在本節中,我們描述一種這樣的方法,似然比檢驗(likelihoodratio test)。似然比檢驗背後的基本思想是計算兩個對比分布下數據的似然性。那麼似然比更高的分布就越適合擬合。或者,可以計算兩個可能性的比率,或者等效地計算比率的對數R,該比率為正或負,取決於更好的分布;如果出現相同的情況,則為零。為了在分布之間做出確定的選擇,我們需要一個對數似然比,該對數似然比必須足夠正或負,以至於不可能由接近零的真實結果引起的機會波動所致。為了定量判斷R的觀測值是否足夠遠離零,我們需要知道預期波動的大小,即我們需要知道R上的標準偏差σ。我們可以使用Vuong提出的方法根據我們的數據進行估算。該方法給出一個p值,該值告訴我們R的觀察符號是否具有統計意義。如果此p值較小(例如p <0.1),則觀察到的信號不太可能是波動的偶然結果,此時該信號能可靠地指示哪個模型更適合數據。另一方面,如果p大,則符號不可靠,並且測試不贊成使用任何一個模型。(In such cases,methods exist which can directly compare two distributions against one anotherand which are considerably easier to implement than the KS test. In thissection we describe one such method, the likelihood ratio test.The basic idea behind thelikelihood ratio test is to compute the likelihood of the data under twocompeting distributions. The one with the higher likelihood is then the betterfit. Alternatively one can calculate the ratio of the two likelihoods, or equivalentlythe logarithm R of the ratio, which is positive or negative depending on whichdistribution is better, or zero in the event of a tie.In order to make a firm choicebetween distributions we need a log likelihood ratio that is sufficientlypositive or negative that it could not plausibly be the result of a chancefluctuation from a true result that is close to zero. To make a quantitativejudgment about whether the observed value of Ris sufficiently far from zero,we need to know the size of the expected fluctuations, i.e., we need to knowthe standard deviation σon R. This we can estimatefrom our data using a method proposed by Vuong. This method gives a p-value that tells us whetherthe observed sign of R is statistically significant. If this p-value is small (say p < 0.1) then it is unlikely thatthe observed sign is a chance result of fluctuations and the sign is a reliableindicator of which model is the better fit to the data. If p is large on the other hand,the sign is not reliable and the test does not favor either model over theother.)
6.實際操作
對於步驟1-3,可以通過powerlaw python package進行powerlaw fit,篩選出較為合適的數據集,對於步驟4可以利用matlab進行p-value檢驗,對於步驟5可以通過powerlaw python package進行檢驗
利用powerlaw python package進行數據擬合
import powerlaw
data = []
fit = powerlaw.Fit(data, discrete=True)
fit.power_law.alpha
fit.power_law.sigma
powerlaw.plot_pdf(data,color='r')
powerlaw.plot_cdf(data,color='g')
powerlaw.plot_ccdf(data,color='b')
powerlaw.plot_pdf(data,linear_bins=True,color='r')
powerlaw.plot_ccdf(data,color='r')
powerlaw.plot_cdf(data,color='r')
fit.plot_pdf()
fit.power_law.plot_pdf()
識別數據中的scaling range
a. 選擇用於擬合的觀測數據集的子集
b. 確定Xmin估計值
c. 利用KS統計量,選擇使得數據子集與擬合模型之間KS距離最小的Xmin
fit.xmin#優化後的Xmin
fit.power_law.D#基於Xmin的數據子集與擬合模型之間KS最小距離
對比不同分布
這裡採用的是前文所述的LRT(loglikelihood ratiostest)
R,p = fit.distribution_compare(『power_law』,『exponential』)
如果一個powerlaw擬合模型不能優於exponential擬合模型,那麼需要非常注意,源數據分布不大可能服從冪律分布。