1
$\begingroup$

I'm trying to evaluate the ks statistic for the gev distribution with ks.test stat function in R. I read the help a few times and remained puzzled as to what the order has to be and can't find any documentation to explain this. The help examples only include a few distributions. This is my short code:

library(fExtremes)
pvec <- runif(2000)
n <- 30 #size of sample
#GEV params
shape <- 0.49; loc <- 36; scale <- 43
#simulate series
gev_sim <- gevSim(model = list(xi = shape, mu = loc, beta = scale), n = n, seed = NULL)
#checking the KS test
tmp <- ks.test(gev_sim,"pgev",loc,scale,shape)
tmp

The P value here for the ks stat is obviously erroneous because the data came from the gev dist. the seemingly best results I get when using this order: loc,scale,shape. But I would like to know what is the correct order.

Also, I would like to know the same for the General Pareto Distribution:

tmp <- ks.test(gp_sim,"pgpd",shape,loc,scale)

Where, loc is the threshold param.

Improve this question
$\endgroup$

1 Answer 1

Reset to default
1
$\begingroup$

According to the documentation of pgev in RStudio's help, the function takes the parameters in the same order as gevSim. To use your nomenclature, the order in which the command is expecting you to input the parameters is pgev(q, xi = shape, mu = loc, beta = scale). Using your code, but changing the order (and also removing the unused pvec), it looks like this:

library(fExtremes)
n <- 30 #size of sample
#GEV params
shape <- 0.49; loc <- 36; scale <- 43
#simulate series
gev_sim <- gevSim(model = list(xi = shape, mu = loc, beta = scale), n = n, seed = NULL)
#checking the KS test
tmp <- ks.test(gev_sim, "pgev", shape, loc, scale)
tmp

and the results I get from one run of the above:

    One-sample Kolmogorov-Smirnov test

data:  gev_sim
D = 0.13785, p-value = 0.5715
alternative hypothesis: two-sided

To further make sure that this is correct, you can increase the sample size (which I think you should anyway, a sample of 30 is unnecessarily small), and check to see whether the resulting data visually seems to fit the corresponding theoretical distribution:

n <- 5000
shape <- 0.49; loc <- 36; scale <- 43
sim_data <- gevSim(model = list(xi = shape, mu = loc, beta = scale), n=n, seed = NULL)
tmp <- ks.test(sim_data, "pgev", shape, loc, scale); tmp

hist(sim_data, freq = FALSE, breaks = 400, col = "steelblue", xlim = c(min(sim_data), quantile(sim_data, 0.995)),
     main = "Comparison of simulated GEV data and theoretical distribution",
     xlab = "x", ylab = "Density")
curve(dgev(x, xi = shape, mu = loc, beta = scale), col = "red", add = TRUE, lwd = 2)
legend("topright", 
       legend=c("Simulated data" , "Density"),  
       lty = c(NA,1), 
       lwd = c(NA,2),
       col = c("black", "red"), 
       border=c("black",NA), 
       fill = c("steelblue",NA), 
       bty = "n", xpd = TRUE, merge = TRUE)

with the output being:

    One-sample Kolmogorov-Smirnov test

data:  sim_data
D = 0.012908, p-value = 0.3754
alternative hypothesis: two-sided

![enter image description here

The case for the GPD is identical to the GEV, the parameters are fed in the function in the exact same way:

gpd_shape <- 0.4; gpd_loc <- 20; gpd_scale <- 27.6
sim_data_gpd <- gpdSim(model = list(xi = gpd_shape, mu = gpd_loc, beta = gpd_scale), n = n, seed = NULL)
tmp_gpd <- ks.test(sim_data_gpd, "pgpd", gpd_shape, gpd_loc, gpd_scale); tmp_gpd
hist(sim_data_gpd, freq = FALSE, breaks = 400, col = "yellow", xlim = c(min(sim_data_gpd), quantile(sim_data_gpd, 0.995)),
     main = "Comparison of simulated GPD data and theoretical distribution",
     xlab = "x", ylab = "Density")
curve(dgpd(x, xi = gpd_shape, mu = gpd_loc, beta = gpd_scale), col = "black", add = TRUE, lwd = 2)
legend("topright", 
       legend=c("Simulated data" , "Density"),  
       lty = c(NA,1), 
       lwd = c(NA,2),
       col = c("black", "black"), 
       border=c("black",NA), 
       fill = c("yellow",NA), 
       bty = "n", xpd = TRUE, merge = TRUE)

resulting in this output:

    One-sample Kolmogorov-Smirnov test

data:  sim_data_gpd
D = 0.010488, p-value = 0.6413
alternative hypothesis: two-sided

enter image description here

Hope this is helpful.

Improve this answer
$\endgroup$
5
  • $\begingroup$ +1 Thank you for a nice illustration of how to confirm what statistical software is doing through testing and simulation. $\endgroup$
    – whuber ♦
    Commented Aug 31, 2018 at 16:18
  • $\begingroup$ Thank you, @whuber! In some runs, the p-value would be at about 0.09, which would reject the null at a sig. level of $\alpha = 0.1$, but I suspect this to be due to the KS test being sensitive to heavy-tailed distributions (but this would require further examination, and I suspect it would only sabotage the purpose of this answer). $\endgroup$
    – Emil
    Commented Aug 31, 2018 at 20:41
  • $\begingroup$ @Emil Perfect, Thanks. I used 30 because my real data is about that and I wanted to examine the P. value sensitivity to those sizes of samples. The detailed answer is very helpful. Cheers. $\endgroup$
    – asher
    Commented Sep 1, 2018 at 8:52
  • $\begingroup$ @asher Ah, my bad then, I thought your data was purely simulated. Regardless, it is useful to keep in mind my comment above about the KS test's sensitivity to heavy tails and its sometimes wrongful rejection of the null (as it did when I was experimenting). In those cases, additional visual inspection (e.g. a QQ plot) or a different test might also be good alternatives. $\endgroup$
    – Emil
    Commented Sep 1, 2018 at 9:07
  • $\begingroup$ @Emil thanks again. For the sake of discussion: I'm not sure how to add a png to the comments but I'm getting the exact opposite impression. I model extreme discharge at rivers and they are all heavy tails. I get KS p.values close to one when obviously, seen from the plot, it is not from the theoretical dist.. Still trying to figure out how to interpret GoF tests... $\endgroup$
    – asher
    Commented Sep 2, 2018 at 8:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.