Goodness-of-fit for Lomax distribution I have some data n > 3000 https://drive.google.com/file/d/1gwB_U_TOX-IQHZJJDX-WeErLzrZZFoXu/view?usp=sharing (Third column) that I believe based on my physical theory should follow a Lomax distribution.
$$
f(x) = \frac{ks^k}{(k + x)^{s+1}} ~~~~~x>0
$$
The parameters of the model are determined by theory and are
$$
s, k = 7.82*10^{-6}, 2.21
$$
The goal is to twofold.
First, I need to determine whether this is a good model for the data particularly near the right tail. Second, I need to compare this model with an exponential fit to the data. To that end, I first used expparetotest from LREP package which test whether data comes from a Pareto distribution and the results confirm that the distribution of data is Pareto and not exponential. I also use the "Jackson" test from the package Renext which gives me the following results.
$statistic
[1] 3.326906

$df
[1] 3185

$p.value
[1] 0

I then use pareto2_test_ad from agop to test if the data is Pareto II, which I have read is the same as Lomax. However, I am getting the following results.
pareto2_test_ad(my_data, 7.8e-6)
W = 467.93, p-value < 2.2e-16

Which if I understand correctly means data is not coming from the specified distribution.
Additionally, I am using result = gpdAd(lowFreq, bootstrap = TRUE, bootnum = 300) from eva
and I get
$statistic
[1] 17.55134

$p.value
[1] 0.003311258

$theta
       Scale        Shape 
7.517058e-07 1.478888e+00 

My questions:

*

*Is my approach correct? Is AD test a good test to determine whether my data comes from Lomax distribution.

*Why do the estimates of the parameters are different from mine?

*Why does p-value change with bootnum?


You could fit the distribution using the gamlss package in R and inspect diagnostic plots. To check whether the fitted parameters are consistent with the theoretical parameters, you could use Wald tests or confidence intervals for the parameters (see below).
Let's fit the distribution and check the diagnostics (I assume the data is loaded and named dat):
library(gamlss)
library(gamlss.dist)
# Fit Lomax distribution
mod <- gamlss(X0~1, family = PARETO2o, data = dat)

plot(mod) # Diagnostic plots
wp(mod, ylim.all = c(1)) # Worm-plot (detrended Q-Q plot)
dtop(mod) # Detrended Owen's plot


Here are the worm plot (upper panel) and Owen's plot (lower panel):

The Q-Q plot and the worm plot as well as Owen's plot show that the Lomax is not a good fit for the data, especially the right tail.
The fitted parameters can be retrieved by:
exp(coef(mod, "mu"))
exp(coef(mod, "sigma"))

So $\hat{\sigma} = 0.6768762$ (this is $k$ in your formula) and $\hat{\mu}= 5.09676\times 10^{-7}$ ($s$ in your formula). To check whether these are consistent with $\sigma = 2.21$ and $\mu = 7.82\times 10^{-6}$, let's do a joint Wald test for both coefficients:
Vmat <- vcov(mod) # Covariance matrix
coefmat <- matrix(c(coef(mod, "mu"), coef(mod, "sigma")), 2, 1) # Coefficient matrix
Lmat <- matrix(c(1, 0, 0, 1), 2, 2) # Hypothesis matrix
cmat <- matrix(c(log(2.21), log(7.82e-6)), 2, 1) # Hypothesized values
# Wald-statistic
Fval <- (1/2)*(t((Lmat%*%coefmat - cmat))%*%solve(Lmat%*%Vmat%*%t(Lmat))%*%(Lmat%*%coefmat - cmat))
# P-value
1 - pf(Fval, 2, df.residual(mod))
[1] 0

The $p$-value is practically $0$ indicating that the model is not consistent with the theoretical values.
Finally, let's calculate the confidence intervals for the parameters:
exp(confint(mod))
                         2.5 %       97.5 %
mu.(Intercept)    4.483221e-07 5.794263e-07
sigma.(Intercept) 6.359687e-01 7.204151e-01

The $95\,\%$ confidence interval for $\mu$ is $(4.48\times 10^{-7};5.79\times 10^{-7})$ and for $\sigma$ it's $(0.636; 0.720)$. The theoretical values are not part of the intervals, reaffirming our conclusion that the fitted model is not compatible with the theoretical paramter values.
