Goodness of fit in R

Question

I need to check for goodness of fit using K-S test in R. I have a dataset (containing 50 data points) and a non-standard continuous probability distribution. Should I peform a one-sample or two-sample test ? In case of one sample test, I don't know if R allows non-standard probability distributions as one of the arguments in ks.test function. If R does allow, I'd like to know how to go about writing the code. If I need to perform a two sample test, then what should be the two datasets that I need to provide as arguments in ks.test function. Also, in both cases, do I have to specify or consider fixed values of the parameters of the distribution ? Is the choice arbitrary ?

Answer 1

According to the manual, ks.test will allow any "continuous (cumulative) distribution function," providing a natural use of a one-sample test in your case for a pre-specified distribution including its parameter values. You just have to write an R function that describes how your "non-standard" distribution function increases from 0 to 1 as a function of the independent variable, and present that function (or a character string naming it) as the argument "y".

If you are estimating parameters of the distribution from the data, however, then the usual distribution of the KS test statistic is no longer valid. As @Glen_b points out, for testing the null hypothesis that the data follow the distribution in such a case, you have to examine the distribution of the KS test statistic (maximum difference between theoretical and empirical distribution functions) under the null hypothesis by generating simulated data sets from the hypothesized distribution. I couldn't find, in a quick look, an R package that accomplishes this in the general case. I did find a recent paper by G.C. Blain that includes examples of R code that accomplishes this task for gamma and normal distributions. I've copied the gamma distribution example below; it should be fairly straightforward to adapt to your distribution.

### Lilliefors test for the 2-parameter gamma distribution
Ns=500000 # number of synthetic samples
n=50 # sample size (given by the user)
shape= 3 # shape parameter (given by the user)
eta=30 # scale parameter (given by the user)
x=matrix(NA,n,1)
lilliefors=matrix(NA,Ns,1)
pos=matrix(1: n, n, 1)/n
for (i in 1:Ns){
    x[,1]=rgamma(n,shape,1/beta)
    A=log(mean(x))-((sum(log(x)))/n)
    alfali=(1/(4*A))*(1+sqrt(1+(4*A/3)))
    etali=mean(x)/alfali
    probpar[,1]=pgamma(sort(x), alfali, 1/betali, lower.tail = TRUE, log.p = FALSE)
    Dmax=max(abs(pos- probpar))
    lilliefors[i,1]=Dmax}
NKScrit5=quantile(lilliefors, probs=0.95) # 5% significance level
NKScrit10=quantile(lilliefors, probs=0.90) # 10% significance level
format(NKScrit5, digits=3)
format(NKScrit10, digits=3)

If you need coding help, Stack Overflow is a better place to look. This site tries to be focused on statistical issues.