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=Bpower,
nbins=c(11, 5), maxProcessor = 1)
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.
CI=FALSE: if TRUE confidence intervals based on the number of MC simulation runs are found.
list.with.everything: this can be a list that has all the functions like pnull, rnull, etc. Usually this was generated by a call to Rgof::case_studies.
B=1000 the number of simulation runs.
It is also possible to pre-define the desired confidence level and run the power routine until this level is achieved, using the routine gof_power_adaptive:
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)/sum(x)/10- 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=Bpower, nbins=c(11, 5), maxProcessor = 1)
#> KS K AD CvM W Wassp1 NS-B NS-A BJ l-P s-P
#> 0.5 0.07 0.04 0.02 0.06 0.06 0.03 0.02 0.04 0.02 0.02 0.02
#> 0.6 0.08 0.03 0.03 0.05 0.05 0.00 0.01 0.01 0.06 0.03 0.03Note 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
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=Bpower, maxProcessor = 1)
#> KS K AD CvM W ZA ZK ZC Wassp1 ES-l-P ES-s-P EP-l-P EP-s-P
#> 0 0.05 0.1 0.09 0.1 0.08 0.08 0.08 0.11 0.12 0.04 0.04 0.03 0.04
#> 1 1.00 1.0 1.00 1.0 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
#> ES-l-L ES-s-L EP-l-L EP-s-L
#> 0 0.04 0.05 0.06 0.04
#> 1 1.00 1.00 1.00 1.00Composite 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=Bpower, maxProcessor=1)
#> KS K AD CvM W ZA ZK ZC Wassp1 ES-l-P ES-s-P EP-l-P EP-s-P
#> 0 0.97 0.99 0.02 0.01 0.03 0.56 0.01 0.00 1.00 0.08 0.08 0.05 0.07
#> 1 0.97 0.99 0.06 0.06 0.08 0.59 0.00 0.03 0.99 0.05 0.03 0.08 0.06
#> ES-l-L ES-s-L EP-l-L EP-s-L
#> 0 0.07 0.08 0.06 0.10
#> 1 0.07 0.06 0.08 0.08ralt = 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=Bpower, maxProcessor=1)
#> KS K AD CvM W ZA ZK ZC Wassp1 ES-l-P ES-s-P EP-l-P EP-s-P
#> 2 1.00 1 0.53 0.54 0.60 0.90 0.19 0.11 0.99 0.19 0.18 0.13 0.18
#> 50 0.95 1 0.02 0.06 0.07 0.51 0.01 0.01 0.55 0.05 0.04 0.04 0.05
#> ES-l-L ES-s-L EP-l-L EP-s-L
#> 2 0.20 0.20 0.13 0.20
#> 50 0.05 0.04 0.03 0.07Power with
Running other tests
It is very easy for a user to add other goodness-of-fit tests to the package.
Example
Say we wish to use tests that are variants 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^p dF(x)\] For continuous data we have the routine
newTScont = function(x, pnull, param) {
Fx=sort(pnull(x))
n=length(x)
out = c(sum(abs( (2*1:n-1)/2/n-Fx )),sum(sqrt(abs( (2*1:n-1)/2/n-Fx ))))
names(out) = c("CvM alt 1", "CvM alt 2")
out
}This routine has to have three or four arguments x, pnull, param and (optionally) TSextra. x is the data and pnull a function that finds the cdf at x. param has to be the estimated parameters in the case of a composite null hypothesis, and is ignored in the case of a simple null hypothesis. Tsextra has to be a list of additional items needed for the calculation of the test statistic, if any.
Note that the return object has to be a named vector.
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)
#> maxProcessor set to 1 because parallel overhead would dominate this short computation
#>
#> Rgof goodness-of-fit tests
#> Data: continuous | n = 500 | B = 5000
#>
#> method statistic p.value
#> CvM alt 1 7.3942 0.3510
#> CvM alt 2 57.3140 0.3262Say 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=3), 3),
Range=c(0,1), B=Bpower, maxProcessor = 1)
#> CvM alt 1 CvM alt 2
#> 0 0.11 0.10
#> 0.25 0.91 0.91
#> 0.5 1.00 1.00for discrete data we will write the routine using Rcpp:
(For reasons to avoid issues with CRAN submission this routine is already part of Rgof)
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
NumericVector newTSdisc(IntegerVector x,
Function pnull,
NumericVector param,
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), Fx(k);
double tmp;
TS.names() = methods;
Fx=pnull(param);
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;
}The routine has to have four or five arguments x, pnull, param, vals and (optionally) TSextra. 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.
As an example we will test whether some data comes from a Binomial distribution with n=10 trials. The success parameter p will be estimated from the data:
vals=0:10
pnull = function(p) pbinom(0:10, 10, p)
rnull = function(p) table(c(0:10,rbinom(10000, 10, p)))-1
phat=function(x) sum(0:10*x)/100000
x = rnull(0.5)
gof_test(x, vals, pnull, rnull, phat=phat, TS=Rgof::newTSdisc)
#> maxProcessor set to 1 because parallel overhead would dominate this short computation
#>
#> Rgof goodness-of-fit tests
#> Data: discrete | n = 10000 | B = 5000
#>
#> method statistic p.value
#> CvM alt 0.001353 0.8000
#> l-P 6.947000 0.6426
#> s-P 6.947000 0.6426
#> l-L 6.938600 0.6435
#> s-L 6.938600 0.6435
#>
#> Estimated parameters:
#> [1] 0.50249For the power calculations we will consider data that is actually a mixture of two Binomials:
ralt = function(tau=0) {
x=rbinom(5000, 10, 0.5-tau)
y=rbinom(5000, 10, 0.5+tau)
table(c(0:10,x,y))-1
}
gof_power(pnull, vals, rnull, ralt,
TS=Rgof::newTSdisc, phat=phat,
param_alt=round(seq(0, 0.05, length=3), 3),
B=Bpower, maxProcessor = 1)
#> CvM alt
#> 0 0.05
#> 0.025 0.23
#> 0.05 1.00If the new routine can also calculate p values use the argument With.p.value=TRUE to indicate that. In that case the return object should be a (vector of) p values.
Adjusted p values for Several Tests
As no single test can be relied upon to consistently have good power, it is reasonable to employ several of them. We could then reject the null hypothesis if any of the tests does so, that is, if the smallest p-value is less than the desired type I error probability \(\alpha\).
This procedure clearly suffers from the problem of simultaneous inference, and the true type I error probability will be much larger than \(\alpha\). It is however possible to adjust the p value so it does achieve the desired \(\alpha\). This can be done as follows:
We generate a number of data sets under the null hypothesis. Generally about 1000 will be sufficient. Then for each simulated data set we apply the tests we wish to include, and record the smallest p value. Here is an example. Say the null hypothesis specifies a uniform \([0.1]\) and a sample size of 250.
pnull=function(x) punif(x)
rnull=function() runif(250)
pvals=matrix(0,1000,16)
for(i in 1:1000)
pvals[i, ]=Rgof::gof_test(rnull(), NA, pnull,
rnull,B=Btest)$p.valuesNote this is not run because of CRAN time constraints. The resulting matrix is in Rgof::pvaluecdf.
Next we find the smallest p value in each run for two selections of four methods. One is the selection found to be best above, namely the methods by Wilson, Anderson-Darling, Zhang’s ZC and a chi square test with a small number of bins and using Pearson’s formula. As a second selection we use the methods by Kolmogorov-Smirnov, Kuiper, Anderson-Darling and Cramer-vonMises. It can be checked that for this null hypothesis these methods are highly correlated.
colnames(pvals)=names(Rgof::gof_test(rnull(), NA, pnull, rnull,B=10)$p.values)
p1=apply(pvals[, c("W", "ZC", "AD", "ES-s-P" )], 1, min)
p2=apply(pvals[, c("KS", "K", "AD", "CvM")], 1, min)Next we find the empirical distribution function for the two sets of p values and draw their graphs. We also add the curve for the cases of four identical tests and the case of four independent tests, which of course is the Bonferroni correction. The data for the cdf is in the inst/extdata directory of the package
tmp=Rgof::pvaluecdf
Tests=factor(c(rep("Identical Tests", nrow(tmp)),
rep("Correlated Selection", nrow(tmp)),
rep("Best Selection", nrow(tmp)),
rep("Independent Tests", nrow(tmp))),
levels=c("Identical Tests", "Correlated Selection",
"Best Selection", "Independent Tests"),
ordered = TRUE)
dta=data.frame(x=c(tmp[,1],tmp[,1],tmp[,1],tmp[,1]),
y=c(tmp[,1],tmp[,3],tmp[,2],1-(1-tmp[,1])^4),
Tests=Tests)
ggplot2::ggplot(data=dta, ggplot2::aes(x=x,y=y,col=Tests))+
ggplot2::geom_line(linewidth=1.2)+
ggplot2::labs(x="p value", y="CDF")+
ggplot2::scale_color_manual(values=c("blue","red", "Orange", "green"))
Here is how to find these adjusted p values with Rgof:
x=rnull()
Rgof::gof_test_adjusted_pvalue(x, NA, pnull, rnull,
B=c(1000,500), maxProcessor = 1)
#> p values of individual tests:
#> KS : 0.444
#> K : 0.14
#> AD : 0.547
#> CvM : 0.59
#> W : 0.25
#> ZA : 0.707
#> ZK : 0.402
#> ZC : 0.622
#> ES-l-P : 0.1164
#> ES-s-P : 0.1274
#> EP-l-P : 0.1994
#> EP-s-P : 0.0048
#> ES-l-L : 0.1001
#> ES-s-L : 0.1351
#> EP-l-L : 0.2174
#> EP-s-L : 0.0067
#> adjusted p value of combined tests: 0.054Weighted 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)-pnorm(-3))/(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)
Rgof::gof_power(pnull, NA, rnull, ralt, w=w, param_alt = c(0,0.2), Range=c(-3,3),B=Bpower, maxProcessor = 1)
#> KS K CvM AD
#> 0 0.06 0.07 0.04 0.03
#> 0.2 1.00 0.98 1.00 1.00It 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 a simulation study to verify that they work in the case under consideration.
Hybrid Tests: Null Distribution Available Only via Simulation
Sometimes it is possible to generate data from the distribution under the null hypothesis, but its cdf can not be calculated. In this case the routines hybrid_test and hybrid_power generate a (usually large) Monte Carlo data set under the null hypothesis and run the two-sample tests of the package R2sample on the data and the Monte Carlo data set. The p values are found with the parametric bootstrap. If parameters are estimated, the Monte Carlo data set is generated with the parameters estimated from the data, and in each bootstrap run the parameters are estimated again from the bootstrap data set, so the dependence between the data and the Monte Carlo data is reproduced. These routines require the package R2sample. If the cdf can be calculated, gof_test should be used, as it generally has higher power. The chi-square tests use the chi-square approximation and treat the Monte Carlo data as an independent sample, so their p values are only approximate when parameters are estimated.
Say we wish to test whether data comes from a normal distribution with mean and standard deviation estimated, pretending that we can not evaluate the normal cdf:
rnull=function(p) rnorm(200, p[1], p[2])
phat=function(x) c(mean(x), sd(x))
x=rt(200, 4)
Rgof::hybrid_test(x, rnull=rnull, phat=phat, nMC=5, B=Bpower, maxProcessor=1)
#>
#> Rgof hybrid goodness-of-fit tests
#> Data: continuous | n = 200 | nMC = 1000 | B = 100
#>
#> method statistic p.value
#> KS 0.0920 0.0000
#> Kuiper 0.1720 0.0000
#> CvM 0.3616 0.0000
#> AD 2.2415 0.0000
#> ZA -3.2780 0.0000
#> ZK 4.6310 0.0000
#> ZC -3.2650 0.0000
#> Wassp1 0.2753 0.0000
#> EppsS 18.8390 0.0000
#> HHG 4.7110 0.0000
#> ES large 31.0800 0.0195
#> ES small 12.6600 0.0130
#> EP large 33.0000 0.1041
#> EP small 21.8400 0.0094
#>
#> Estimated parameters:
#> [1] -0.07757627 1.62349243The argument nMC is the size of the Monte Carlo data set, or, if it is at most 10, its size as a multiple of the sample size.
Power Graphs
It is easy to generate power graphs with the routine plot_power. Say we wish to generate a power graph for the case of a uniform distribution if the data actually comes from a linear one:
rnull=function() runif(n)
pnull=function(x) punif(x)
TSextra(qnull=function(x) qunif(x))
ralt=function(a=0.3) {
if(a==0) return(runif(n))
(-(1-a)+sqrt((1-a)^2+4*a*runif(n)))/2/a
}
pwr=Rgof::gof_power(pnull, NA, rnull, ralt,
param_alt =seq(0.01, 0.5, length=21),
Range=c(0,1))
Case Studies
The package includes the routine run.studies, which provides 110 case studies each for the continuous and the discrete case. This allows the user to easily compare the power of the methods included in the package to a different one of their choice.
The 120 studies fall into 5 types. In the first four, with 20 cases each, we have data sets that differ in exactly one population parameter: either mean, variance, skewness or kurtosis. In type 5 (=“mixed”) any of the four parameters can be different. The first 20 of these have fixed parameters and the last 20 have parameters estimated from the data.
As an example consider a type 1 (=“mean”) study. Study #1 has one data set from a uniform distribution on \([0,1]\). Therefore we have mean=0.5, variance=1/12, skewness=0 and kurtosis=1.8. Next we find a distribution with the same variance, skewness and kurtosis, but with a mean of 0.53. Their densities look as follows:
f=Rgof::case_studies("uniform", "mean", n=1e4)
#> Error:
#> ! 'case_studies' is not an exported object from 'namespace:Rgof'
x=0:100/100
dta=data.frame(x=c(x,x),
y=c(f$dnull(x), Rgof::funs_list[["cont"]][["uniform"]][["mean"]]$density(x)),
Distribution=rep(c("Uniform",
"Alternative"), each=101))
#> Error:
#> ! object 'f' not found
ggplot(dta, aes(x, y, col=Distribution)) +
geom_line() +
scale_y_continuous(limits = c(0, 2)) +
labs(y="Density")
#> Error in `ggplot()`:
#> ! could not find function "ggplot"These alternative distributions are designed so that for a true type I error of \(\alpha=0.05\) the best method as a power of between \(70\%\) and \(90\%\).
The package includes two example “new” tests, both simple chi-square tests, for either continuous (myTS_cont) or discrete data (myTS_disc). Either can return the value of the test statistic, or a p value, depending on what is specified in TSextra.
These cases can be run easily with the routine run.studies:
TSextra=list(nbins=5, statistic=FALSE)
Rgof::run.studies(1:2, 1:2,
TS=Rgof::myTS_cont, TSextra=TSextra,
With.p.value = TRUE)
#> Average rank of a method:
#> ZC ZA ZK AD EP-l-P ES-l-P EP-l-L EP-s-P ES-l-L ES-s-P EP-s-L
#> 100.0 94.4 88.9 79.2 61.8 61.1 56.2 55.6 54.9 47.9 46.5
#> ES-s-L Chi_C5 Wassp1 K CvM W KS
#> 41.7 34.7 33.3 30.6 22.9 22.9 17.4
#> Chi_C5 KS K AD CvM W ZA ZK ZC Wassp1
#> uniform/mean 0.614 0.448 0.203 0.739 0.529 0.118 0.825 0.813 0.840 0.560
#> linear/mean 0.298 0.182 0.096 0.508 0.253 0.091 0.842 0.788 0.885 0.249
#> uniform/mean 0.109 0.082 0.148 0.310 0.090 0.130 0.802 0.757 0.873 0.122
#> linear/mean 0.099 0.047 0.104 0.167 0.047 0.100 0.699 0.641 0.793 0.065
#> ES-l-P ES-s-P EP-l-P EP-s-P ES-l-L ES-s-L EP-l-L EP-s-L
#> uniform/mean 0.431 0.684 0.431 0.684 0.426 0.682 0.426 0.682
#> linear/mean 0.470 0.287 0.423 0.334 0.368 0.288 0.382 0.319
#> uniform/mean 0.130 0.120 0.165 0.120 0.138 0.114 0.169 0.114
#> linear/mean 0.239 0.113 0.178 0.117 0.197 0.110 0.154 0.113Arguments of run.studies:
- studies: either the name(s) or position(s) of the desired studies in
the list of all studies
- types: one of 5 types of departures from the null.
- data_type=“cont” for continuous data or “disc” for discrete data.
- TS: new test statistic
- TSextra: a list passed to TS
- With.p.value =FALSE set to TRUE if the new method finds p value(s).
TS should then return only those
- n: sample size
- alpha=0.05: true type I error probability
- B=1000: number of simulation runs
Arguments of routine to calculate new test:
Continuous data
- x: the data set - pnull: routine to calculate the cdf
- p: vector with parameter estimates, for composite hypotheses.
- TSextra: (optional): a list passed to the routine with any object
needed to do the calculation.
Discrete data Same as for continuous data plus vals, a vector of the values of the discrete random variable.
The output of the routine has to be a named vector with either the test statistic(s) or the p value(s).
The case studies are as follows. As discussed above, there are 20 studies with one distribution chosen from the usual list, and the second so it matches the three of the four parameters mean, variance, skewness and kurtosis but differs in the fourth. The distributions are:
- uniform
- linear
- quadratic
- uniform-betabump
- sine
- betaaa
- beta2a
- triangular
- uniformmixture
- betamixture
- normal-pure
- normal-mix1
- normal-mix2
- normal-mix3
- normal-mix4
- gamma1
- gamma2
- gammamixture
- exponential-mixture
- noncentral
For the exact definitions run case.studies:
Rgof::case_studies("noncentral")[2:4]
#> Error:
#> ! 'case_studies' is not an exported object from 'namespace:Rgof'In addition there are another 30 case studies where most of the parameters are different:
uniform.linear\(\hspace{3cm}\) U[0,1] vs a linear model on [0,1] with slope s.
uniform.quadratic\(\hspace{2.5cm}\) U[0,1] vs a quadratic model with vertex at 0.5 and some curvature a.
uniform.bump\(\hspace{3.1cm}\) U[0,1] vs U[0,1]+N(0.5,0.05).
uniform.sine\(\hspace{3.3cm}\) U[0,1] vs U[0,1]+Sine wave
beta22.betaaa\(\hspace{3cm}\) Beta(2,2) vs Beta(a,a)
beta22.beta2a\(\hspace{3cm}\) Beta(2,2) vs Beta(2,a)
uniform.uniformmixture\(\hspace{1.9cm}\) U[0,1] vs. \(\alpha\)U[0,1/2]+(1-\(\alpha\))U[1/2,1]
uniform.betamixture\(\hspace{2.5cm}\) U[0,1] vs. \(\alpha\)U[0,1/2]+(1-\(\alpha\))Beta(2,2)
uniform.triangular\(\hspace{2.9cm}\) U[0,1] vs. triangular
double.exponential\(\hspace{3.5cm}\) with different rates.
normal.onestretch\(\hspace{3.1cm}\) N(0,1) vs 1/2N(0,1)+1/2N(0, \(\sigma\))
normal.t\(\hspace{4.2cm}\) N(0,1) vs t(df)
normal.outlier1\(\hspace{3.1cm}\) N(0,1) vs N(0,1)+N(5, 1)
normal.outlier2\(\hspace{3.1cm}\) N(0,1) vs N(0,1)+U[-5,-3]+U[3,5]
normal.normalmixture\(\hspace{3.1cm}\) N(0,1) vs N(-a,1)+N(a,1)
exponential.gamma\(\hspace{2.3cm}\) Exp(1) vs Gamma(1,b)
exponential.weibull\(\hspace{2.5cm}\) Exp(1) vs Weibull(1,b)
exponential.bump\(\hspace{2.7cm}\) Exp(1) vs Exp(1)+N(0.5,0.05)
chisquare.noncentral\(\hspace{2.3cm}\) \(\chi^2(5)\) vs. \(\chi^2(5, \tau)\)
gamma.gammamixture\(\hspace{2.3cm}\) \(\Gamma(2, 3)\) vs \((1-\alpha)\Gamma(2, 3)+ \alpha\Gamma(3,3)\)
The last 10 case studies include parameter estimation:
normal.t.est\(\hspace{3.1cm}\) \(N(\mu, \sigma)\) vs \(t(5)\), with mean and standard deviation estimated.
exponential.weibull.est\(\hspace{3.1cm}\) \(Exp(\lambda)\) vs Weibull(1.2, 1), rate estimated
trunc.exponential.linear.est\(\hspace{3.1cm}\) 0-1 truncate exponential with rate \(\lambda\) vs linear, rate estimated
exponential.gamma.est\(\hspace{3.1cm}\) \(Exp(\lambda)\) vs \(\Gamma(1.2, 1)\), rate estimated
normal.cauchy.est\(\hspace{3.1cm}\) \(N(0,\sigma)\) vs Cauchy, both truncated to the interval \([-3.5, 3.5\), variance estimated
betaa1.betaab.est\(\hspace{3.1cm}\) Beta(a, 1) vs Beta(1, 1.2), shape1 estimated
betaa1.betaaa.est\(\hspace{3.1cm}\) Beta(a, 1) vs Beta(1.2, 1.2), shape1 estimated
trans-beta.est\(\hspace{3.1cm}\) \(g(x)=(\exp(x)-1)/(e-1);g(X)~Beta(a,1)\) vs. \(h(x)=(\exp(hx)-1)/(\exp(a)-1);h(X)~Beta(1,1)\), a estimated
double.exponential.est\(\hspace{3.1cm}\) double exponential rate \(\lambda\) vs. double gamma(1.6, 1), rate estimated
pareto.est\(\hspace{3.1cm}\) \(f(x)=\alpha/(x+1)^{\alpha + 1}\) vs Gamma(x, 1.07, 8), \(\alpha\) estimated.
betamixture1.est\(\hspace{3.1cm}\) \(f(x)=a*dbeta(x, 1, 1)+(1-a)*dbeta(x, 2, 2)\), mixing ratio a estimated.
“betamixture2.est\(\hspace{3.1cm}\) \(f(x)=a*dbeta(x, 3, 1)+(1-a)*dbeta(x, 1, 3)\), mixing ratio a estimated.
betamixture3.est\(\hspace{3.1cm}\) \(f(x)=a*dbeta(x, 2, 2)+(1-a)*dbeta(x, 2, 5)\), mixing ratio a estimated.
normalmixture.est\(\hspace{3.1cm}\) f(x)=a*pnorm(x)+(1-a)*pnorm(x, 2), mixing ratio a estimated.
laplace1.est\(\hspace{3.1cm}\) f(x, p)=function(x, p=c(0, 1)) exp(-abs(x - p[1]) / p[2]) / (2 * p[2]), location and scale estimated. Alternative is a mixture of normals.
laplace2.est\(\hspace{3.1cm}\) f(x, p)=exp(-abs(x - p[1]) / p[2]) / (2 * p[2]), alternative is a double gamma distribution.
truncexp1.est\(\hspace{3.1cm}\) f(x,l,a)=dexp(x, l)/pexp(a, l), upper limit a and rate l estimated.
truncexp2.est\(\hspace{3.1cm}\) f(x, l)=dexp(x, l)/(pexp(2, l)-pexp(1, l)), rate l estimated.
gammamixture1.est\(\hspace{3.1cm}\) f(x, a)=a*dgamma(x, 1, 1)+(1-a)*dgamma(x, 3, 1), mixing ratio a estimated.
gammamixture2.est\(\hspace{3.1cm}\) f(x, a)=a*dgamma(x, 3, 1)+(1-a)*dgamma(x, 10, 1), mixing ratio a estimated.
For the discrete versions of case studies 55-57 the parameters are estimated with the grouped-data maximum likelihood estimator, which maximizes \(\sum_j n_j\log\left(F(b_j)-F(b_{j-1})\right)\) over the bins \((b_{j-1},b_j]\) (for 57 the upper limit is the upper edge of the last nonempty bin). The continuous estimators applied to the bin midpoints (the median, the maximum) take only a few distinct values, and the parametric bootstrap then does not hold its level.
Generally a user will wish to apply their test to all the case studies. This can be done easily with:
run.studies can also be used to run the studies for the included methods but for different values of the sample size n and/or alpha. Note, however, that these studies are designed so that for the pre-chosen sample sizes and \(\alpha=0.05\) the power of the best method is generally larger than \(70\%\) and less than \(90\%\), so as to make comparisons meaningfull.
or a user can recreate the file Rgof::power_study: