poweRlaw (version 0.20.2)

compare_distributions: Vuong's test for non-nested models

Description

Since it is possible to fit power law models to any data set, it is recommended that different distributions are considered. A standard technique is to use Vuong's test, which a likelihood ratio test for model selection using the Kullback-Leibler criteria. The test statistic, R, is the ratio of the log-likelihoods of the data between the two competing models. The sign of R indicates which model is better. Since the value of R is obviously subject to error, we use the method proposed by Vuong, 1989.

Usage

compare_distributions(d1, d2)

Value

This function returns [object Object],[object Object],[object Object],[object Object]

Details

This function compares two models. The null hypothesis is that both classes of distributions are equally far from the true distribution. If this is true, the log-likelihood ratio should (asymptotically) have a Normal distribution with mean zero. The test statistic is the sample average of the log likelihood ratio, standardized by a consistent estimate of its standard deviation. If the null hypothesis is false, and one class of distributions is closer to the "truth", this test statistic goes to +/-infinity with probability 1, indicating the better-fitting class of distributions.

References

Vuong, Quang H. (1989): "Likelihood Ratio Tests for Model Selection and Non-Nested Hypotheses", Econometrica 57: 307--333.

Examples

########################################################
# Example data                                         #
########################################################
x = 3:10; xmin = 2

########################################################
##Continuous power law                                 #
########################################################
m1 = conpl$new(x)
m1$setXmin(xmin)
est1 = estimate_pars(m1)
m1$setPars(est1$pars)

########################################################
##Exponential                                          #
########################################################
m2 = conexp$new(x)
m2$setXmin(xmin)
est2 = estimate_pars(m2)
m2$setPars(est2$pars)

########################################################
##Vuong's test                                         #
########################################################
comp = compare_distributions(m1, m2)
plot(comp)