poweRlaw (version 0.20.3)

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 alternative distributions are considered. A standard technique is to use Vuong's test. This is 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 esimated, we use the method proposed by Vuong, 1989 to select the model. 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", the test statistic goes to +/-infinity with probability 1, indicating the better-fitting class of distributions.

Usage

compare_distributions(d1, d2)

Value

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

References

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

Examples

########################################################
# Example data                                         #
########################################################
x = rpldis(100, xmin=2, alpha=3)

########################################################
##Continuous power law                                 #
########################################################
m1 = conpl$new(x)
m1$setXmin(estimate_xmin(m1))

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

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