Nothing
R Package Rgof"
In Rgof: 1d Goodness of Fit Tests
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
error = TRUE
)
library(Rgof)
Bsim = c(100, 200) #Number of Simulation Runs
The package Rgof brings together a number of routines for the goodness-of-fit problem for univariate data. We have a data set $\pmb{x}$, and we want to test whether it was generated by the probability distribution F.
The highlights of this package are:
- it runs a large number of different tests simultaneously.
- it can run several tests and find an adjusted p value for the combination of all the tests.
- it works for continuous and for discrete (aka histogram) data.
- it allows for parameter estimation.
- it uses Rcpp and parallel programming to be very fast.
- it includes routines that make it very easy to estimate the power of the different tests for specific situations.
- it allows for the possibility that the sample size was a draw from a Poisson random variable.
- it allows the user to add more tests.
- it includes a routine to draw the usual power graph.
set.seed(123)
Note all runs of the test routine are done with B=1000 and all runs of the power routines with arguments B=c(500, 500), maxProcessor = 2 in order to pass devtools::check().
The Methods
-
Continuous Data
-
Kolmogorov Smirnov (KS) [@massey1951], [@kolmogorov1933], [@smirnov1948]
- Kuiper (K) [@kuiper1960]
- Anderson-Darling (AD) [@anderson1952], [@anderson1954]
- Cramer-vonMises (CvM) [@anderson1962]
- Wilson (W)
- Zhang's methods (ZA, ZK, ZC) [@zhang2002]
- Wasserstein p=1 (Wassp1) [@wasserstein1999]
For all of these tests the distribution of the test statistic under the null hypothesis is found via simulation.
- Eight variations of chi square tests, using the formulas by Pearson or based on the likelihood ratio test, bins of equal size or equal probability and both a large number and a small number of bins, 50 and 10 by default. The p values are found using the usual chi square approximation to the null distribution. If parameters are estimated this is done either via the method of minimum chi square [@berkson1980] or via a user-provided estimator. In all cases bins are combined until all of them have an expected count of at least 5.
There is a very large literature on chi square tests, the oldest of the goodness of fit tests. For a survey see [@rolke2020].
- Discrete (Histogram) Data
All the methods above are also implemented for discrete data, except for Zhang's tests, which have no discrete analog.
It is worth noting that these discrete versions are based on the theoretical ideas of the tests and not on the actual formula of calculation for the continuous case. The test statistics can therefore be different even when applied to the same data. For example, the Anderson-Darling test is based on the distance measure
$$A^2=n\int_{-\infty}^{\infty} \frac{(\hat{F}(x)-F(x))^2}{F(x)(1-F(x))}dF(x) $$
where $F$ is the theoretical distribution function under the null hypothesis and $\hat{F}$ is the empirical distribution function. In the case of continuous data it can be shown that
$$A^2=-n-\frac1n\sum_{i=1}^n (2i-1)\left(\log F(x_i) +\log[1-F(x_{n+1-i})\right)$$
However, for discrete data we have
$$A^2=n\sum_{i=1}^k \frac{(\hat{F}(x_i)-F(x_i))^2}{F(x_i)(1-F(x_i))}\left(F(x_i)-F(x_{i-1}\right)$$
with $F(x_0)=0$.
In the continuous case $\hat{F}$ is a step function but $F$ is continuous, and therefore $A^2>0$. In the discrete case however$A^2=0$ is possible. This shows that the two cases are fundamentally different.
As for continuous data null distributions are found using simulation. In fact in the case of discrete data none of the tests has a known distribution for the test statistic under the null hypothesis.
- Four variations of chi square tests, using the formulas by Pearson and log-likelihood as well as a large number and a small number of bins. Again the routine combines bins until all have expected counts greater than 5, and the chi square approximation is used to find p values. The combination of bins is done in such a way that the bins remain of equal size as much as possible.
These methods can be used for both discrete and histogram data. The main difference between these two is that discrete data has (a countable) number of possible values whereas histogram data has possible ranges of values (the bins). The only method directly affected by this difference is Wassp1, which requires actual values. All other methods ignore the vals argument.
Testing
Discrete (Histogram) Data/Model
Simple Null Hypothesis
We generate a data set of size 1000 from a Binomial distribution with n=20 and success probability 0.5, and then test $H_0:F=Bin(20, 0.5)$.
vals=0:20 #possible values of random variable
pnull=function() pbinom(0:20, 20, 0.5) # cumulative distribution function (cdf)
rnull = function() table(c(0:20, rbinom(1000, 20, 0.5)))-1
# generate data under the null hypothesis, make sure that vector of counts has same length as vals, possibly 0.
- Null Hypothesis is true
x = rnull()
# Basic Test
gof_test(x, vals, pnull, rnull, B=1000)
#Test with adjusted overall p value
gof_test_adjusted_pvalue(x, vals, pnull, rnull, B=c(1000, 500))
- Null Hypothesis is false
x = table(c(0:20, rbinom(1000, 20, 0.55)))-1
#true p is 0.55, not 0.5
# Basic Test
gof_test(x, vals, pnull, rnull, B=1000, doMethod = "all")$p.value
#Test with adjusted overall p value
gof_test_adjusted_pvalue(x, vals, pnull, rnull, B=c(1000, 500))
Arguments of gof_test for discrete data/model:
-
x: vector with counts (histogram heights). Should have a number for each value of vals, possibly 0.
-
vals: all possible values of discrete random variable, that is all x with $P(X=x)>0$
-
pnull: function to find values of cumulative distribution function for each value of vals. Function has no arguments.
-
rnull: function to generate data from true density. Function has no arguments. Function needs to insure that output is a vector with same length as vals.
-
B=5000: number of simulation runs
-
w: function to find importance sampling weights, if needed
-
phat: function to estimate parameters
-
TS: function to find values of user-supplied test statistics
-
TSextra: a list that is passed to TS if any additional info is required.
-
nbins=c(50, 10): number of bins for chi square tests. The first one is already given by the data in the discrete case, the for the second bins are joined.
-
rate=0, if not 0 sample size is assumed to have come from a Poisson random variable with rate "rate".
-
minexpcount=5, minimal expected counts for chi square tests.
-
ChiUsePhat=TRUE, if TRUE uses phat for parameter estimation. If false uses method of minimum chi square
-
maxProcessors=1 if greater than 1 number of cores for parallel processing.
-
doMethods="all" names of methods to include
The arguments of gof_test_adjusted_pvalue for discrete data/model are the same, except that the number of simulation runs B is two numbers. The first is used for estimating the individual p values, the second for the adjustment.
Random Sample Size
In some fields like high energy physics it is common that the sample size is not fixed but a random variable drawn from a Poisson distribution with a known rate. Our package easily accomodates that:
rnull = function() table(c(0:20, rbinom(rpois(1, 650), 20, 0.5)))-1
x = rnull()
gof_test(x, vals, pnull, rnull, rate=650, B=1000)$p.value
Composite Null Hypothesis
We generate a data set of size 1000 from a binomial distribution with n=20 and success probability p, and then test F=Bin(20, .). p is estimated from data.
vals=0:20
pnull=function(p=0.5) pbinom(0:20, 20, ifelse(p>0&&p<1, p, 0.5))
rnull = function(p=0.5) table(c(0:20, rbinom(1000, 20, p)))-1
phat = function(x) sum(0:20*x)/sum(x)/20
- Null Hypothesis is true
x = table(c(0:20, rbinom(1000, 20, 0.5)))-1
gof_test(x, vals, pnull, rnull, phat=phat, B=1000)$p.value
- Null Hypothesis is true
x = table(c(0:20, rbinom(1000, 20, 0.55)))-1
# p is not 0.5, but data is still from a binomial distribution with n=20
gof_test(x, vals, pnull, rnull, phat=phat, B=1000)$p.value
- Null Hypothesis is false
x = table(c(rep(0:20, 5), rbinom(1000-21*5, 20, 0.53)))
# data has to many small and large values to be from a binomial
gof_test(x, vals, pnull, rnull, phat=phat, B=1000)$p.value
The estimation of the parameter(s) in the case of the chi square tests is done either by using the function phat or via the minimum chi square method. The routine uses a general function minimizer. If there are values of the parameter that are not possible this can lead to warnings. It is best to put a check into the pnull function to avoid this issue. As an example the function pnull above checks that the success probability p is in the interval $(0,1)$.
Histogram Data
A variant of discrete data sometimes encountered is data given in the form of a histogram, that is as a set of bins and their counts. The main distinction is that discrete data has specific values, for example the non-negative integers for a Poisson distribution, whereas histogram data has ranges of numbers, the bins. It turns out that, though, that the only method that requires actual values is Wassp1, and for that method one can use the midpoint of the intervals.
As an example consider the following case: we have histogram data and we want to test whether it comes from an exponential rate 1 distribution, truncated to the interval 0-2:
rnull = function() {
y = rexp(2500, 1) # Exp(1) data
y = y[y<2][1:1500] # 1500 events on 0-2
bins = 0:40/20 # binning
hist(y, bins, plot=FALSE)$counts # find bin counts
}
x = rnull()
bins = 0:40/20
vals = (bins[-1]+bins[-21])/2
pnull = function() {
bins = 1:40/20
pexp(bins, 1)/pexp(2, 1)
}
gof_test(x, vals, pnull, rnull)$p.value
Continuous Data
Simple Hypothesis
pnull = function(x) pnorm(x)
rnull = function() rnorm(1000)
TSextra = list(qnull=function(x) qnorm(x)) #optional quantile function used by chi square tests and Wassp1 test.
- Null Hypothesis is true
x = rnorm(1000)
#Basic Tests
gof_test(x, NA, pnull, rnull, B=1000, TSextra=TSextra)$p.value
#Adjusted p value
gof_test_adjusted_pvalue(x, NA, pnull, rnull, B=c(1000,500), TSextra=TSextra)
- Null Hypothesis is false
x = rnorm(1000, 0.5)
gof_test(x, NA, pnull, rnull, B=1000, TSextra=TSextra)$p.value
Composite Hypothesis - One Parameter
pnull = function(x, p=0) pnorm(x, p)
TSextra = list(qnull = function(x, p=0) qnorm(x, p))
rnull = function(p) rnorm(1000, p)
phat = function(x) mean(x)
- Null Hypothesis is true
x = rnorm(1000)
gof_test(x, NA, pnull, rnull, phat=phat, TSextra=TSextra, B=1000)$p.value
- Null Hypothesis is true
x = rnorm(1000, 0.5)
gof_test(x, NA, pnull, rnull, phat=phat, TSextra=TSextra)$p.value
- Null Hypothesis is false
x = rnorm(1000, 0.5, 2)
gof_test(x, NA, pnull, rnull, phat=phat, TSextra=TSextra, B=1000)$p.value
Composite Hypothesis - Multiple Parameters
pnull = function(x, p=c(0, 1)) pnorm(x, p[1], ifelse(p[2]>0, p[2], 0.001))
TSextra = list(qnull = function(x, p=c(0, 1)) qnorm(x, p[1], ifelse(p[2]>0, p[2], 0.001)))
rnull = function(p=c(0, 1)) rnorm(1000, p[1], ifelse(p[2]>0, p[2], 0.001))
phat = function(x) c(mean(x), sd(x))
- Null Hypothesis is true
x = rnorm(1000)
gof_test(x, NA, pnull, rnull, phat=phat, TSextra=TSextra, B=1000)$p.value
- Null Hypothesis is true
x = rnorm(1000, 0.5)
gof_test(x, NA, pnull, rnull, phat=phat, TSextra=TSextra, B=1000)$p.value
- Null Hypothesis is true
x = rnorm(1000, 0.5, 2)
gof_test(x, NA, pnull, rnull, phat=phat, TSextra=TSextra, B=1000)$p.value
- Null Hypothesis is false
x = rt(1000, 2)
gof_test(x, NA, pnull, rnull, phat=phat, TSextra=TSextra, B=1000)$p.value
Power Estimation
For estimating the power of the various tests one also has to provide the routine ralt, which generates data under the alternative hypothesis:
Discrete Data/Model
Simple Null Hypothesis
vals = 0:10
pnull = function() pbinom(0:10, 10, 0.5)
rnull =function () table(c(0:10, rbinom(100, 10, 0.5)))-1
ralt =function (p=0.5) table(c(0:10, rbinom(100, 10, p)))-1
P=gof_power(pnull, vals, rnull, ralt,
param_alt=seq(0.5, 0.6, 0.02), B=Bsim, nbins=c(11, 5))
plot_power(P, "p", Smooth=FALSE)
In all cases the arguments are the same as for gof_test. In addition we now have
-
ralt: a routine with one parameter that generates data under some alternative hypothesis.
-
param_alt: values to be passed to ralt. This allows the calculation of the power for many different values.
-
alpha=0.05 type I error probability for tests.
-
B=c(1000, 1000) the first number is the number of simulation runs for power estimation and the second the number of runs to be used to find the null distribution.
Composite Null Hypothesis
vals = 0:10
pnull = function(p=0.5) pbinom(0:10, 10, ifelse(0<p&p<1,p,0.001))
rnull = function (p=0.5) table(c(0:10, rbinom(100, 10, ifelse(0<p&p<1,p,0.001))))-1
phat = function(x) sum(0:10*x)/1000
- Null Hypothesis is true
ralt =function (p=0.5) table(c(0:10, rbinom(100, 10, p)))-1
gof_power(pnull, vals, rnull, ralt, c(0.5, 0.6), phat=phat,
B=Bsim, nbins=c(11, 5), maxProcessors = 2)
Note that power estimation in the case of a composite hypothesis (aka with parameters estimated) is much slower than the simple hypothesis case.
- Null Hypothesis is false
ralt =function (p=0.5) table(c(rep(0:10, 2), rbinom(100, 10, p)))
gof_power(pnull, vals, rnull, ralt, 0.5, phat=phat,
B=Bsim, nbins=c(11, 5), maxProcessors = 2)
Continuous Data/Model
Simple Null Hypothesis
pnull = function(x) pnorm(x)
TSextra = list(qnull = function(x) qnorm(x))
rnull = function() rnorm(100)
ralt = function(mu=0) rnorm(100, mu)
gof_power(pnull, NA, rnull, ralt, c(0, 1), TSextra=TSextra, B=Bsim)
Composite Null Hypothesis
pnull = function(x, p=c(0,1)) pnorm(x, p[1], ifelse(p[2]>0, p[2], 0.01))
TSextra = list(qnull = function(x, p=c(0,1)) qnorm(x, p[1], ifelse(p[2]>0, p[2], 0.01)))
rnull = function(p=c(0,1)) rnorm(500, p[1], p[2])
ralt = function(mu=0) rnorm(100, mu)
phat = function(x) c(mean(x), sd(x))
gof_power(pnull, NA, rnull, ralt, c(0, 1), phat= phat,
TSextra=TSextra, B=Bsim, maxProcessor=2)
ralt = function(df=1) {
# t distribution truncated at +- 5
x=rt(1000, df)
x=x[abs(x)<5]
x[1:100]
}
gof_power(pnull, NA, rnull, ralt, c(2, 50), phat=phat,
Range=c(-5,5), TSextra=TSextra, B=Bsim, maxProcessor=2)
Running other tests
Its is very easy for a user to add other goodness-of-fit tests to the package. This can be done by editing the routines TS_cont and/or TS_disc, which are located in the folder inst/examples in the Rgof library folder. Or a user can write their own version of these files.
Example
Say we wish to use a test that is a variant of the Cramer-vonMises test, using the integrated absolute difference of the empirical and the theoretical distribution function:
$$\int_{-\infty}^{\infty} \vert F(x) - \hat{F}(x) \vert dF(x)$$
For continuous data we have the routine
newTScont = function(x, Fx) {
Fx=sort(Fx)
n=length(x)
out = sum(abs( (2*1:n-1)/2/n-Fx ))
names(out) = "CvM alt"
out
}
This routine has to have two arguments x and Fx. Note that the return object has to be a named vector. The object TSextra can be used to provide further information to the TS routine, if necessary.
Then we can run this test with
pnull = function(x) punif(x)
rnull = function() runif(500)
x = rnull()
Rgof::gof_test(x, NA, pnull, rnull, TS=newTScont)
Say we want to find the power of this test when the true distribution is a linear:
ralt = function(slope=0) {
if(slope==0) y=runif(500)
else y=(slope-1+sqrt((1-slope)^2+4*slope* runif(500)))/2/slope
}
gof_power(pnull, NA, rnull, ralt, TS=newTScont, param_alt=round(seq(0, 0.5, length=5), 3), Range=c(0,1), B=Bsim)
for discrete data we write will the routine using Rcpp:
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
NumericVector newTSdisc(IntegerVector x,
NumericVector Fx,
NumericVector vals) {
Rcpp::CharacterVector methods=CharacterVector::create("CvM alt");
int const nummethods=methods.size();
int k=x.size(), n, i;
NumericVector TS(nummethods), ecdf(k);
double tmp;
TS.names() = methods;
n=0;
for(i=0;i<k;++i) n = n + x[i];
ecdf(0) = double(x(0))/double(n);
for(i=1;i<k;++i) {
ecdf(i) = ecdf(i-1) + x(i)/double(n);
}
tmp = std::abs(ecdf[0]-Fx(0))*Fx(0);
for(i=1;i<k;++i)
tmp = tmp + std::abs(ecdf(i)-Fx(i))*(Fx(i)-Fx(i-1));
TS(0) = tmp;
return TS;
}
Again the routine has to have three arguments x, Fx and vals and the output vector has to have names.
Note that one drawback of writing the routine in Rcpp is that it is then not possible to use multiple processors.
vals=1:50/51
pnull = function() (1:50)/50
rnull = function() c(rmultinom(1, 500, rep(1/50,50)))
x = rnull()
gof_test(x, vals, pnull, rnull, TS=newTSdisc)
and for power calculations we can run
ralt = function(slope=0) {
if(slope==0) p=rep(1/50, 50)
else p=diff(slope * (0:50/50)^2 + (1 - slope) * 0:50/50)
c(rmultinom(1, 500, p))
}
gof_power(pnull, vals, rnull, ralt, TS=newTSdisc, param_alt=round(seq(0, 0.5, length=5), 3), B=Bsim)
Weighted Data
Sometimes the data/model uses importance sampling weights. This can be done as follows. Say we want to test whether the data comes from a standard normal distribution, truncated to [-3,3] and with weights from a t distribution with 3 degrees of freedom:
$H_0: F=N(0,1)$, $X\sim t(3)$
df=3
pnull=function(x) pnorm(x)/(2*pnorm(3)-1)
rnull=function() {x=rt(2000, df);x=x[abs(x)<3];sort(x[1:1000])}
w=function(x) (dnorm(x)/(2*pnorm(3)-1))/(dt(x,df)/(2*pt(3,df)-1))
x=sort(rnull())
plot(x, w(x), type="l", ylim=c(0, 2*max(w(x))))
ralt=function(m=0) {x=rt(2000,df)+m;x=x[abs(x)<3];sort(x[1:1000])}
set.seed(111)
gof_power(pnull, NA, rnull, ralt, w=w, param_alt = c(0,0.2), Range=c(-3,3),B=Bsim)
It should be noted that these tests are quite sensitive to the size of the weights and to the sample size, so one should always do simulation study to verify that they work in the case under consideration.
References
Try the Rgof package in your browser
Any scripts or data that you put into this service are public.
Rgof documentation built on Oct. 9, 2024, 1:06 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set( collapse = TRUE, comment = "#>", error = TRUE )
library(Rgof) Bsim = c(100, 200) #Number of Simulation Runs
The package Rgof brings together a number of routines for the goodness-of-fit problem for univariate data. We have a data set $\pmb{x}$, and we want to test whether it was generated by the probability distribution F.
The highlights of this package are:
- it runs a large number of different tests simultaneously.
- it can run several tests and find an adjusted p value for the combination of all the tests.
- it works for continuous and for discrete (aka histogram) data.
- it allows for parameter estimation.
- it uses Rcpp and parallel programming to be very fast.
- it includes routines that make it very easy to estimate the power of the different tests for specific situations.
- it allows for the possibility that the sample size was a draw from a Poisson random variable.
- it allows the user to add more tests.
- it includes a routine to draw the usual power graph.
set.seed(123)
Note all runs of the test routine are done with B=1000 and all runs of the power routines with arguments B=c(500, 500), maxProcessor = 2 in order to pass devtools::check().
The Methods
-
Continuous Data
-
Kolmogorov Smirnov (KS) [@massey1951], [@kolmogorov1933], [@smirnov1948]
- Kuiper (K) [@kuiper1960]
- Anderson-Darling (AD) [@anderson1952], [@anderson1954]
- Cramer-vonMises (CvM) [@anderson1962]
- Wilson (W)
- Zhang's methods (ZA, ZK, ZC) [@zhang2002]
- Wasserstein p=1 (Wassp1) [@wasserstein1999]
For all of these tests the distribution of the test statistic under the null hypothesis is found via simulation.
- Eight variations of chi square tests, using the formulas by Pearson or based on the likelihood ratio test, bins of equal size or equal probability and both a large number and a small number of bins, 50 and 10 by default. The p values are found using the usual chi square approximation to the null distribution. If parameters are estimated this is done either via the method of minimum chi square [@berkson1980] or via a user-provided estimator. In all cases bins are combined until all of them have an expected count of at least 5.
There is a very large literature on chi square tests, the oldest of the goodness of fit tests. For a survey see [@rolke2020].
- Discrete (Histogram) Data
All the methods above are also implemented for discrete data, except for Zhang's tests, which have no discrete analog.
It is worth noting that these discrete versions are based on the theoretical ideas of the tests and not on the actual formula of calculation for the continuous case. The test statistics can therefore be different even when applied to the same data. For example, the Anderson-Darling test is based on the distance measure
$$A^2=n\int_{-\infty}^{\infty} \frac{(\hat{F}(x)-F(x))^2}{F(x)(1-F(x))}dF(x) $$ where $F$ is the theoretical distribution function under the null hypothesis and $\hat{F}$ is the empirical distribution function. In the case of continuous data it can be shown that
$$A^2=-n-\frac1n\sum_{i=1}^n (2i-1)\left(\log F(x_i) +\log[1-F(x_{n+1-i})\right)$$ However, for discrete data we have
$$A^2=n\sum_{i=1}^k \frac{(\hat{F}(x_i)-F(x_i))^2}{F(x_i)(1-F(x_i))}\left(F(x_i)-F(x_{i-1}\right)$$
with $F(x_0)=0$.
In the continuous case $\hat{F}$ is a step function but $F$ is continuous, and therefore $A^2>0$. In the discrete case however$A^2=0$ is possible. This shows that the two cases are fundamentally different.
As for continuous data null distributions are found using simulation. In fact in the case of discrete data none of the tests has a known distribution for the test statistic under the null hypothesis.
- Four variations of chi square tests, using the formulas by Pearson and log-likelihood as well as a large number and a small number of bins. Again the routine combines bins until all have expected counts greater than 5, and the chi square approximation is used to find p values. The combination of bins is done in such a way that the bins remain of equal size as much as possible.
These methods can be used for both discrete and histogram data. The main difference between these two is that discrete data has (a countable) number of possible values whereas histogram data has possible ranges of values (the bins). The only method directly affected by this difference is Wassp1, which requires actual values. All other methods ignore the vals argument.
Testing
Discrete (Histogram) Data/Model
Simple Null Hypothesis
We generate a data set of size 1000 from a Binomial distribution with n=20 and success probability 0.5, and then test $H_0:F=Bin(20, 0.5)$.
vals=0:20 #possible values of random variable pnull=function() pbinom(0:20, 20, 0.5) # cumulative distribution function (cdf) rnull = function() table(c(0:20, rbinom(1000, 20, 0.5)))-1 # generate data under the null hypothesis, make sure that vector of counts has same length as vals, possibly 0.
- Null Hypothesis is true
x = rnull() # Basic Test gof_test(x, vals, pnull, rnull, B=1000) #Test with adjusted overall p value gof_test_adjusted_pvalue(x, vals, pnull, rnull, B=c(1000, 500))
- Null Hypothesis is false
x = table(c(0:20, rbinom(1000, 20, 0.55)))-1 #true p is 0.55, not 0.5 # Basic Test gof_test(x, vals, pnull, rnull, B=1000, doMethod = "all")$p.value #Test with adjusted overall p value gof_test_adjusted_pvalue(x, vals, pnull, rnull, B=c(1000, 500))
Arguments of gof_test for discrete data/model:
-
x: vector with counts (histogram heights). Should have a number for each value of vals, possibly 0.
-
vals: all possible values of discrete random variable, that is all x with $P(X=x)>0$
-
pnull: function to find values of cumulative distribution function for each value of vals. Function has no arguments.
-
rnull: function to generate data from true density. Function has no arguments. Function needs to insure that output is a vector with same length as vals.
-
B=5000: number of simulation runs
-
w: function to find importance sampling weights, if needed
-
phat: function to estimate parameters
-
TS: function to find values of user-supplied test statistics
-
TSextra: a list that is passed to TS if any additional info is required.
-
nbins=c(50, 10): number of bins for chi square tests. The first one is already given by the data in the discrete case, the for the second bins are joined.
-
rate=0, if not 0 sample size is assumed to have come from a Poisson random variable with rate "rate".
-
minexpcount=5, minimal expected counts for chi square tests.
-
ChiUsePhat=TRUE, if TRUE uses phat for parameter estimation. If false uses method of minimum chi square
-
maxProcessors=1 if greater than 1 number of cores for parallel processing.
-
doMethods="all" names of methods to include
The arguments of gof_test_adjusted_pvalue for discrete data/model are the same, except that the number of simulation runs B is two numbers. The first is used for estimating the individual p values, the second for the adjustment.
Random Sample Size
In some fields like high energy physics it is common that the sample size is not fixed but a random variable drawn from a Poisson distribution with a known rate. Our package easily accomodates that:
rnull = function() table(c(0:20, rbinom(rpois(1, 650), 20, 0.5)))-1 x = rnull() gof_test(x, vals, pnull, rnull, rate=650, B=1000)$p.value
Composite Null Hypothesis
We generate a data set of size 1000 from a binomial distribution with n=20 and success probability p, and then test F=Bin(20, .). p is estimated from data.
vals=0:20 pnull=function(p=0.5) pbinom(0:20, 20, ifelse(p>0&&p<1, p, 0.5)) rnull = function(p=0.5) table(c(0:20, rbinom(1000, 20, p)))-1 phat = function(x) sum(0:20*x)/sum(x)/20
- Null Hypothesis is true
x = table(c(0:20, rbinom(1000, 20, 0.5)))-1 gof_test(x, vals, pnull, rnull, phat=phat, B=1000)$p.value
- Null Hypothesis is true
x = table(c(0:20, rbinom(1000, 20, 0.55)))-1 # p is not 0.5, but data is still from a binomial distribution with n=20 gof_test(x, vals, pnull, rnull, phat=phat, B=1000)$p.value
- Null Hypothesis is false
x = table(c(rep(0:20, 5), rbinom(1000-21*5, 20, 0.53))) # data has to many small and large values to be from a binomial gof_test(x, vals, pnull, rnull, phat=phat, B=1000)$p.value
The estimation of the parameter(s) in the case of the chi square tests is done either by using the function phat or via the minimum chi square method. The routine uses a general function minimizer. If there are values of the parameter that are not possible this can lead to warnings. It is best to put a check into the pnull function to avoid this issue. As an example the function pnull above checks that the success probability p is in the interval $(0,1)$.
Histogram Data
A variant of discrete data sometimes encountered is data given in the form of a histogram, that is as a set of bins and their counts. The main distinction is that discrete data has specific values, for example the non-negative integers for a Poisson distribution, whereas histogram data has ranges of numbers, the bins. It turns out that, though, that the only method that requires actual values is Wassp1, and for that method one can use the midpoint of the intervals.
As an example consider the following case: we have histogram data and we want to test whether it comes from an exponential rate 1 distribution, truncated to the interval 0-2:
rnull = function() { y = rexp(2500, 1) # Exp(1) data y = y[y<2][1:1500] # 1500 events on 0-2 bins = 0:40/20 # binning hist(y, bins, plot=FALSE)$counts # find bin counts } x = rnull() bins = 0:40/20 vals = (bins[-1]+bins[-21])/2 pnull = function() { bins = 1:40/20 pexp(bins, 1)/pexp(2, 1) } gof_test(x, vals, pnull, rnull)$p.value
Continuous Data
Simple Hypothesis
pnull = function(x) pnorm(x) rnull = function() rnorm(1000) TSextra = list(qnull=function(x) qnorm(x)) #optional quantile function used by chi square tests and Wassp1 test.
- Null Hypothesis is true
x = rnorm(1000) #Basic Tests gof_test(x, NA, pnull, rnull, B=1000, TSextra=TSextra)$p.value #Adjusted p value gof_test_adjusted_pvalue(x, NA, pnull, rnull, B=c(1000,500), TSextra=TSextra)
- Null Hypothesis is false
x = rnorm(1000, 0.5) gof_test(x, NA, pnull, rnull, B=1000, TSextra=TSextra)$p.value
Composite Hypothesis - One Parameter
pnull = function(x, p=0) pnorm(x, p) TSextra = list(qnull = function(x, p=0) qnorm(x, p)) rnull = function(p) rnorm(1000, p) phat = function(x) mean(x)
- Null Hypothesis is true
x = rnorm(1000) gof_test(x, NA, pnull, rnull, phat=phat, TSextra=TSextra, B=1000)$p.value
- Null Hypothesis is true
x = rnorm(1000, 0.5) gof_test(x, NA, pnull, rnull, phat=phat, TSextra=TSextra)$p.value
- Null Hypothesis is false
x = rnorm(1000, 0.5, 2) gof_test(x, NA, pnull, rnull, phat=phat, TSextra=TSextra, B=1000)$p.value
Composite Hypothesis - Multiple Parameters
pnull = function(x, p=c(0, 1)) pnorm(x, p[1], ifelse(p[2]>0, p[2], 0.001)) TSextra = list(qnull = function(x, p=c(0, 1)) qnorm(x, p[1], ifelse(p[2]>0, p[2], 0.001))) rnull = function(p=c(0, 1)) rnorm(1000, p[1], ifelse(p[2]>0, p[2], 0.001)) phat = function(x) c(mean(x), sd(x))
- Null Hypothesis is true
x = rnorm(1000) gof_test(x, NA, pnull, rnull, phat=phat, TSextra=TSextra, B=1000)$p.value
- Null Hypothesis is true
x = rnorm(1000, 0.5) gof_test(x, NA, pnull, rnull, phat=phat, TSextra=TSextra, B=1000)$p.value
- Null Hypothesis is true
x = rnorm(1000, 0.5, 2) gof_test(x, NA, pnull, rnull, phat=phat, TSextra=TSextra, B=1000)$p.value
- Null Hypothesis is false
x = rt(1000, 2) gof_test(x, NA, pnull, rnull, phat=phat, TSextra=TSextra, B=1000)$p.value
Power Estimation
For estimating the power of the various tests one also has to provide the routine ralt, which generates data under the alternative hypothesis:
Discrete Data/Model
Simple Null Hypothesis
vals = 0:10 pnull = function() pbinom(0:10, 10, 0.5) rnull =function () table(c(0:10, rbinom(100, 10, 0.5)))-1 ralt =function (p=0.5) table(c(0:10, rbinom(100, 10, p)))-1 P=gof_power(pnull, vals, rnull, ralt, param_alt=seq(0.5, 0.6, 0.02), B=Bsim, nbins=c(11, 5)) plot_power(P, "p", Smooth=FALSE)
In all cases the arguments are the same as for gof_test. In addition we now have
-
ralt: a routine with one parameter that generates data under some alternative hypothesis.
-
param_alt: values to be passed to ralt. This allows the calculation of the power for many different values.
-
alpha=0.05 type I error probability for tests.
-
B=c(1000, 1000) the first number is the number of simulation runs for power estimation and the second the number of runs to be used to find the null distribution.
Composite Null Hypothesis
vals = 0:10 pnull = function(p=0.5) pbinom(0:10, 10, ifelse(0<p&p<1,p,0.001)) rnull = function (p=0.5) table(c(0:10, rbinom(100, 10, ifelse(0<p&p<1,p,0.001))))-1 phat = function(x) sum(0:10*x)/1000
- Null Hypothesis is true
ralt =function (p=0.5) table(c(0:10, rbinom(100, 10, p)))-1 gof_power(pnull, vals, rnull, ralt, c(0.5, 0.6), phat=phat, B=Bsim, nbins=c(11, 5), maxProcessors = 2)
Note that power estimation in the case of a composite hypothesis (aka with parameters estimated) is much slower than the simple hypothesis case.
- Null Hypothesis is false
ralt =function (p=0.5) table(c(rep(0:10, 2), rbinom(100, 10, p))) gof_power(pnull, vals, rnull, ralt, 0.5, phat=phat, B=Bsim, nbins=c(11, 5), maxProcessors = 2)
Continuous Data/Model
Simple Null Hypothesis
pnull = function(x) pnorm(x) TSextra = list(qnull = function(x) qnorm(x)) rnull = function() rnorm(100) ralt = function(mu=0) rnorm(100, mu) gof_power(pnull, NA, rnull, ralt, c(0, 1), TSextra=TSextra, B=Bsim)
Composite Null Hypothesis
pnull = function(x, p=c(0,1)) pnorm(x, p[1], ifelse(p[2]>0, p[2], 0.01)) TSextra = list(qnull = function(x, p=c(0,1)) qnorm(x, p[1], ifelse(p[2]>0, p[2], 0.01))) rnull = function(p=c(0,1)) rnorm(500, p[1], p[2]) ralt = function(mu=0) rnorm(100, mu) phat = function(x) c(mean(x), sd(x)) gof_power(pnull, NA, rnull, ralt, c(0, 1), phat= phat, TSextra=TSextra, B=Bsim, maxProcessor=2)
ralt = function(df=1) { # t distribution truncated at +- 5 x=rt(1000, df) x=x[abs(x)<5] x[1:100] } gof_power(pnull, NA, rnull, ralt, c(2, 50), phat=phat, Range=c(-5,5), TSextra=TSextra, B=Bsim, maxProcessor=2)
Running other tests
Its is very easy for a user to add other goodness-of-fit tests to the package. This can be done by editing the routines TS_cont and/or TS_disc, which are located in the folder inst/examples in the Rgof library folder. Or a user can write their own version of these files.
Example
Say we wish to use a test that is a variant of the Cramer-vonMises test, using the integrated absolute difference of the empirical and the theoretical distribution function:
$$\int_{-\infty}^{\infty} \vert F(x) - \hat{F}(x) \vert dF(x)$$ For continuous data we have the routine
newTScont = function(x, Fx) { Fx=sort(Fx) n=length(x) out = sum(abs( (2*1:n-1)/2/n-Fx )) names(out) = "CvM alt" out }
This routine has to have two arguments x and Fx. Note that the return object has to be a named vector. The object TSextra can be used to provide further information to the TS routine, if necessary.
Then we can run this test with
pnull = function(x) punif(x) rnull = function() runif(500) x = rnull() Rgof::gof_test(x, NA, pnull, rnull, TS=newTScont)
Say we want to find the power of this test when the true distribution is a linear:
ralt = function(slope=0) { if(slope==0) y=runif(500) else y=(slope-1+sqrt((1-slope)^2+4*slope* runif(500)))/2/slope }
gof_power(pnull, NA, rnull, ralt, TS=newTScont, param_alt=round(seq(0, 0.5, length=5), 3), Range=c(0,1), B=Bsim)
for discrete data we write will the routine using Rcpp:
#include <Rcpp.h> using namespace Rcpp; // [[Rcpp::export]] NumericVector newTSdisc(IntegerVector x, NumericVector Fx, NumericVector vals) { Rcpp::CharacterVector methods=CharacterVector::create("CvM alt"); int const nummethods=methods.size(); int k=x.size(), n, i; NumericVector TS(nummethods), ecdf(k); double tmp; TS.names() = methods; n=0; for(i=0;i<k;++i) n = n + x[i]; ecdf(0) = double(x(0))/double(n); for(i=1;i<k;++i) { ecdf(i) = ecdf(i-1) + x(i)/double(n); } tmp = std::abs(ecdf[0]-Fx(0))*Fx(0); for(i=1;i<k;++i) tmp = tmp + std::abs(ecdf(i)-Fx(i))*(Fx(i)-Fx(i-1)); TS(0) = tmp; return TS; }
Again the routine has to have three arguments x, Fx and vals and the output vector has to have names.
Note that one drawback of writing the routine in Rcpp is that it is then not possible to use multiple processors.
vals=1:50/51 pnull = function() (1:50)/50 rnull = function() c(rmultinom(1, 500, rep(1/50,50))) x = rnull() gof_test(x, vals, pnull, rnull, TS=newTSdisc)
and for power calculations we can run
ralt = function(slope=0) { if(slope==0) p=rep(1/50, 50) else p=diff(slope * (0:50/50)^2 + (1 - slope) * 0:50/50) c(rmultinom(1, 500, p)) } gof_power(pnull, vals, rnull, ralt, TS=newTSdisc, param_alt=round(seq(0, 0.5, length=5), 3), B=Bsim)
Weighted Data
Sometimes the data/model uses importance sampling weights. This can be done as follows. Say we want to test whether the data comes from a standard normal distribution, truncated to [-3,3] and with weights from a t distribution with 3 degrees of freedom:
$H_0: F=N(0,1)$, $X\sim t(3)$
df=3 pnull=function(x) pnorm(x)/(2*pnorm(3)-1) rnull=function() {x=rt(2000, df);x=x[abs(x)<3];sort(x[1:1000])} w=function(x) (dnorm(x)/(2*pnorm(3)-1))/(dt(x,df)/(2*pt(3,df)-1)) x=sort(rnull()) plot(x, w(x), type="l", ylim=c(0, 2*max(w(x)))) ralt=function(m=0) {x=rt(2000,df)+m;x=x[abs(x)<3];sort(x[1:1000])} set.seed(111) gof_power(pnull, NA, rnull, ralt, w=w, param_alt = c(0,0.2), Range=c(-3,3),B=Bsim)
It should be noted that these tests are quite sensitive to the size of the weights and to the sample size, so one should always do simulation study to verify that they work in the case under consideration.
References
Try the Rgof package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.