knitr::opts_chunk$set(
echo = TRUE,
error = TRUE,
comment = "##>"
)
seed0=1234
set.seed(seed0)
B=250
maxP=1 #necessary for full reproducibility
library(MDgof) # Ver 1.2.1
library(MD2sample) # Ver 1.4.2##
## Attaching package: 'MD2sample'## The following objects are masked from 'package:MDgof':
##
## case.studies, run.studieswhich.to.run=c(TRUE, TRUE, TRUE, TRUE,
TRUE, TRUE, TRUE, TRUE,
TRUE, TRUE, TRUE,
TRUE, TRUE)
names(which.to.run) =c("GCT", "GDT", "TCT", "TDT",
"GCP", "GDP", "TCP", "TDP",
"HCP", "run1", "run2",
"CompTime1", "CompTime2")pnullC = function (x) {
f = function(x) mvtnorm::pmvnorm(rep(-Inf, length(x)), x)
if (!is.matrix(x))
return(f(x))
apply(x, 1, f)
}
rnullC = function() mvtnorm::rmvnorm(100, c(0, 0))
set.seed(seed0)
x=rnullC()
gof_test(x, pnullC, rnullC, B=B, maxProcessor=maxP, seed=seed0)##>
##> MDgof multivariate goodness-of-fit tests
##> Data: continuous | n = 100 | dimension = 2 | B = 250
##>
##> method statistic p.value
##> qK 0.1740 0.1920
##> qCvM 0.2940 0.0440
##> qAD 2.1220 0.2320
##> FF 0.4610 0.7240
##> RK 0.0495 0.8480
##> EP 18.1000 0.1538dnullC = function (x) {
f = function(x) mvtnorm::dmvnorm(x)
if (!is.matrix(x))
return(f(x))
apply(x, 1, f)
}
gof_test(x, pnullC, rnullC, dnull=dnullC,
B=B, maxProcessor=maxP, seed=seed0)##>
##> MDgof multivariate goodness-of-fit tests
##> Data: continuous | n = 100 | dimension = 2 | B = 250
##>
##> method statistic p.value
##> qK 0.17400 0.1920
##> qCvM 0.29400 0.0440
##> qAD 2.12200 0.2320
##> FF 0.46100 0.7240
##> RK 0.04950 0.8480
##> EP 18.10000 0.1538
##> BB 0.00647 0.8360options(width=65)
MDgof::hybrid_test(x, rnullC, nMC=2, B=B, maxProcessor=maxP, seed=seed0)##> Data is assumed to be continuous##>
##> MD2sample two-sample tests
##> Data: continuous | n.x = 100 | n.y = 200 | dimension = 2 | B = 250
##>
##> method statistic p.value
##> KS 0.13000 0.372
##> K 0.22500 0.132
##> CvM 0.17710 0.232
##> AD 1.11680 0.292
##> NN1 1.18000 0.012
##> NN5 2.30000 0.000
##> AZ 0.06003 0.368
##> BF 0.25570 0.452
##> BG 0.07833 0.728
##> MMD 0.01627 0.340set.seed(seed0)
pnullD=function(x, p) {
f=function(x) sum(dbinom(0:x[1], 10, p[1])*pbinom(x[2], 5, 0.5+0:x[1]/100))
if(!is.matrix(x)) x=rbind(x)
apply(x, 1, f)
}
rnullD=function(p=0.5) {
x=stats::rbinom(1000, 10, p)
y=stats::rbinom(1000, 5, 0.5+x/100)
MDgof::sq2rec(table(x, y))
}
phatD=function(x) {
tx=tapply(x[,3], x[,1], sum)
mean(rep(as.numeric(names(tx)), times=tx))/10
}
set.seed(seed0)
x=rnullD(0.57)
gof_test(x, pnullD, rnullD, phat=phatD,
B=B, maxProcessor=maxP, seed=seed0)##>
##> MDgof multivariate goodness-of-fit tests
##> Data: discrete | n = 1000 | dimension = 2 | B = 250
##>
##> method statistic p.value
##> K 0.04810 0.104
##> AD 0.00127 0.116
##> KL 0.04260 0.016
##> P 89.90000 0.072
##>
##> Estimated parameters:
##> [1] 0.5652x=stats::rbinom(1000, 2, 0.5)
y=stats::rbinom(1000, 2, 0.7+x/100)
head(cbind(x,y))##> x y
##> [1,] 1 2
##> [2,] 0 1
##> [3,] 2 2
##> [4,] 0 1
##> [5,] 0 1
##> [6,] 1 1MDgof::sq2rec(table(x, y))##> [,1] [,2] [,3]
##> [1,] 0 0 30
##> [2,] 1 0 38
##> [3,] 2 0 13
##> [4,] 0 1 107
##> [5,] 1 1 216
##> [6,] 2 1 96
##> [7,] 0 2 121
##> [8,] 1 2 244
##> [9,] 2 2 135set.seed(seed0)
options(width = 65)
x=matrix(runif(500), ncol=5)
y=matrix(runif(500), ncol=5)
MD2sample::twosample_test(x, y, B=B,
maxProcessor=maxP, seed=seed0)##> Data is assumed to be continuous##>
##> MD2sample two-sample tests
##> Data: continuous | n.x = 100 | n.y = 100 | dimension = 5 | B = 250
##>
##> method statistic p.value
##> KS 0.10000 0.6040
##> K 0.19000 0.3120
##> CvM 0.03440 0.3280
##> AD 0.98140 0.5200
##> NN1 0.89000 0.8960
##> NN5 5.11000 0.2680
##> AZ 0.09747 0.0400
##> BF 0.62240 0.0320
##> BG 2.92460 0.0320
##> MMD 0.02420 0.0600
##> FR 0.56830 0.7151
##> NN0 0.22000 0.8644
##> CF1 0.56830 0.7151
##> CF2 1.48400 0.4762
##> CF3 -0.56800 0.7151
##> CF4 1.07750 0.4051
##> Ball 0.01260 0.0236raltD=function(K=0.5) {
x=stats::rbinom(1000, 10, 0.5)
y=stats::rbinom(1000, 5, K+x/100)
MDgof::sq2rec(table(x, y))
}set.seed(seed0)
x = table(sample(1:4, size=1000, replace = TRUE))
y = table(sample(1:4, size=1000, replace = TRUE,
prob=c(1,1.5,1,1)))
vals_x=rep(1:2,2)
vals_y=rep(1:2, each=2)
twosample_test(x, y, vals_x, vals_y,
B=B, seed=seed0)##> Data is assumed to be discrete##> maxProcessor set to 1 because parallel overhead would dominate this short computation##>
##> MD2sample two-sample tests
##> Data: discrete | n.x = 1000 | n.y = 1000 | dimension = 2 | B = 250
##>
##> method statistic p.value
##> KS 0.08400 0e+00
##> K 0.13200 0e+00
##> CvM 1.41150 0e+00
##> AD 5.95920 0e+00
##> NN 10.07200 0e+00
##> AZ 0.03334 0e+00
##> BF 0.08924 0e+00
##> ChiSquare 19.76900 2e-04set.seed(seed0)
raltC = function(s) {
mvtnorm::rmvnorm(100, c(0, 0), sigma=matrix(c(1, s, s, 1), 2, 2))
}
MDgof::gof_power(pnullC, rnullC, raltC, c(0, 0.5), B=B, maxProcessor=maxP, seed=seed0)##> qKS qK qCvM qAD BR FF RK ES EP
##> 0 0.052 0.080 0.032 0.048 0.084 0.084 0.068 0.072 0.084
##> 0.5 0.404 0.208 0.576 0.624 0.088 0.252 0.148 0.700 0.816TSextra=list(Continuous=TRUE,
WithEstimation=FALSE,
Withpvalue=TRUE,
Range=matrix(c(-Inf, Inf, -Inf, Inf), 2, 2))
MDgof::gof_power(pnullC, rnullC, raltC, c(0, 0.5),
TS=newTS, TSextra=TSextra,
With.p.value = TRUE,
B=B, maxProcessor=maxP, seed=seed0)##> ES33 EP33
##> 0 0.056 0.048
##> 0.5 0.596 0.860options(width=65)
MDgof::hybrid_power(rnullC, raltC, c(0, 0.5), B=B, maxProcessor=maxP, seed=seed0)##> KS K CvM AD NN1 NN5 AZ BF BG MMD
##> 0 0.088 0.060 0.076 0.076 0.072 0.052 0.048 0.056 0.052 0.060
##> 0.5 0.116 0.136 0.196 0.276 0.248 0.260 0.220 0.188 0.040 0.284
##> FR NN0 CF1 CF2 CF3 CF4 Ball ES EP
##> 0 0.032 0.092 0.032 0.044 0.032 0.036 0.048 0.056 0.064
##> 0.5 0.164 0.148 0.164 0.104 0.164 0.116 0.080 0.416 0.316set.seed(seed0)
options(width=65)
MDgof::gof_power(pnullD, rnullD, raltD, c(0.5, 0.53),
phat=phatD, B=B, maxProcessor=maxP, seed=seed0)##> KS K CvM AD TV KL H P
##> 0.5 0.024 0.012 0.020 0.024 0.068 0.068 0.060 0.124
##> 0.53 0.928 0.844 0.928 0.936 0.536 0.556 0.468 0.372set.seed(seed0)
f = function (theta = 0, d = 2) {
x = matrix(stats::runif(500), ncol = 5)
if (theta == 0)
return(list(x=x, y=matrix(runif(500), ncol=5)))
cop = copula::frankCopula(theta, 5)
list(x = x, y = copula::rCopula(100, cop))
}
MD2sample::twosample_power(f, c(0, 3.5), B=B,
maxProcessor=maxP, seed=seed0)##> KS K CvM AD NN1 NN5 AZ BF BG MMD
##> 0 0.052 0.056 0.040 0.036 0.056 0.048 0.052 0.052 0.056 0.044
##> 3.5 0.872 0.596 0.988 0.988 0.912 0.992 1.000 0.988 0.060 1.000
##> FR NN0 CF1 CF2 CF3 CF4 Ball
##> 0 0.044 0.052 0.044 0.032 0.044 0.052 0.032
##> 3.5 0.936 0.544 0.936 0.924 0.936 0.916 0.948set.seed(seed0)
MDgof::case.studies(1:3, Continuous = TRUE,
WithEstimation = TRUE,
Dim=5,
ReturnCaseNames = TRUE)[1:10]##> [1] "uniform.uniform-diagonal-n" "uniform.uniform-diagonal-b"
##> [3] "normal-ind.normal-corA" "normal-ind.normal-corB"
##> [5] "t-ind.t-cor" "uniform.Frank"
##> [7] "uniform.Clayton" "uniform.Gumbel"
##> [9] "uniform.Joe" "mix.uniform.Clayton"tp=proc.time()
options(width=70)
TSextra=list(Continuous=TRUE,
WithEstimation=FALSE,
Withpvalue=TRUE)
res=MDgof::run.studies(TS=MDgof::newTS, TSextra=TSextra,
With.p.value = TRUE, B=B,
maxProcessor=maxP, seed=seed0)##> Running case study uniform.uniform-diagonal-n for continuous 2 dimensional data and without parameter estimation##> Running case study uniform.uniform-diagonal-b for continuous 2 dimensional data and without parameter estimation##> Running case study normal-ind.normal-cor for continuous 2 dimensional data and without parameter estimation##> Running case study t-ind.t-cor for continuous 2 dimensional data and without parameter estimation##> Running case study uniform.Frank for continuous 2 dimensional data and without parameter estimation##> Running case study uniform.Clayton for continuous 2 dimensional data and without parameter estimation##> Running case study uniform.Gumbel for continuous 2 dimensional data and without parameter estimation##> Running case study uniform.Galambos for continuous 2 dimensional data and without parameter estimation##> Running case study uniform.HuslerReiss for continuous 2 dimensional data and without parameter estimation##> Running case study uniform.Joe for continuous 2 dimensional data and without parameter estimation##> Running case study mix.Clayton.Gumbel for continuous 2 dimensional data and without parameter estimation##> Running case study mix.uniform.Frank for continuous 2 dimensional data and without parameter estimation##> Running case study mix.Clayton.Frank for continuous 2 dimensional data and without parameter estimation##> Running case study mix.Frank.Gumbel for continuous 2 dimensional data and without parameter estimation##> Running case study mix.Frank.Joe for continuous 2 dimensional data and without parameter estimation##> Running case study normal.normal-shift-one.marginal for continuous 2 dimensional data and without parameter estimation##> Running case study normal.normal-shift-both.marginal for continuous 2 dimensional data and without parameter estimation##> Running case study normal.normal-stretch-one.marginal for continuous 2 dimensional data and without parameter estimation##> Running case study normal.normal-stretch-both.marginal for continuous 2 dimensional data and without parameter estimation##> Running case study normal.stretch-shift-one.marginal for continuous 2 dimensional data and without parameter estimation##> Running case study normal.stretch-shift-two.marginal for continuous 2 dimensional data and without parameter estimation##> Running case study uniform.rotated-uniform.marginal for continuous 2 dimensional data and without parameter estimation##> Running case study uniform.beta-one.marginal for continuous 2 dimensional data and without parameter estimation##> Running case study uniform.beta-two.marginal for continuous 2 dimensional data and without parameter estimation##> Running case study uni-exp-1.uni-exp-l.marginal for continuous 2 dimensional data and without parameter estimation##> Running case study exp-exp-1.exp-exp-l.marginal for continuous 2 dimensional data and without parameter estimation##> Running case study beta-nor-1.beta-nor-mean.marginal for continuous 2 dimensional data and without parameter estimation##> Running case study beta-nor-1.beta-nor-sd.marginal for continuous 2 dimensional data and without parameter estimation##> Running case study beta-beta-2.beta-beta-a.marginal for continuous 2 dimensional data and without parameter estimation##> Running case study beta05.normal.marginal for continuous 2 dimensional data and without parameter estimation##> Average rank of each test across studies:##> ES33 EP33 qKS qK qCvM qAD BB BR KSD FF RK ES EP
##> 10.0 9.3 6.7 6.3 7.0 6.8 2.0 5.1 5.2 8.2 5.4 9.6 9.3proc.time()-tp##> user system elapsed
##> 195.63 4.54 208.32res[1:5, ]##> ES33 EP33 qKS qK qCvM qAD BB
##> uniform.uniform-diagonal-n 0.964 0.948 0.396 0.294 0.566 0.691 0.075
##> uniform.uniform-diagonal-b 0.916 0.928 0.352 0.250 0.525 0.776 0.138
##> normal-ind.normal-cor 0.784 0.944 0.370 0.199 0.584 0.699 0.142
##> t-ind.t-cor 0.480 0.900 0.424 0.225 0.628 0.673 0.108
##> uniform.Frank 0.968 0.948 0.449 0.319 0.616 0.670 0.039
##> BR KSD FF RK ES EP
##> uniform.uniform-diagonal-n 0.498 0.508 0.621 0.411 0.953 0.959
##> uniform.uniform-diagonal-b 0.390 0.379 0.485 0.212 0.946 0.958
##> normal-ind.normal-cor 0.272 0.441 0.629 0.061 0.864 0.861
##> t-ind.t-cor 0.080 0.094 0.564 0.055 0.666 0.838
##> uniform.Frank 0.555 0.384 0.615 0.071 0.835 0.844set.seed(seed0)
tmp=MDgof::case.studies(1, nsample=1000)
set.seed(seed0)
x=tmp$rnull()
tp=proc.time()
a=MDgof::gof_test(x, tmp$pnull, tmp$rnull,
dnull=tmp$dnull, B=B,
maxProcessor=maxP, seed=seed0)
proc.time()-tp##> user system elapsed
##> 41.50 0.56 43.09x=matrix(runif(500), ncol=5)
y=matrix(runif(500), ncol=5)
tp=proc.time()
a=MD2sample::twosample_test(x, y, B=B,
maxProcessor=maxP, seed=seed0)##> Data is assumed to be continuousproc.time()-tp##> user system elapsed
##> 3.50 0.22 3.80pnull = function (x, p) {
f = function(x)
mvtnorm::pmvnorm(rep(-Inf, length(x)), x,
mean=p[1:4],
sigma=matrix(p[5:20], nrow=4))
if (!is.matrix(x)) return(f(x))
apply(x, 1, f)
}
dnull = function (x, p) {
f = function(x)
mvtnorm::dmvnorm(x,
mean=p[1:4],
sigma=matrix(p[5:20], nrow=4))
if (!is.matrix(x)) return(f(x))
apply(x, 1, f)
}
rnull = function(p) mvtnorm::rmvnorm(50,
mean=p[1:4],
sigma=matrix(p[5:20], nrow=4))
phat=function(x) c(apply(x, 2, mean), cov(x))setosa=iris3[,,1]
p <- GGally::ggpairs(
setosa,
progress = FALSE
) +
ggplot2::theme_bw(base_size = 13) +
ggplot2::theme(
axis.text = ggplot2::element_text(size = 10),
strip.text = ggplot2::element_text(size = 12),
strip.text.x = ggplot2::element_text(size = 12),
strip.text.y = ggplot2::element_text(size = 12)
)
MDgof::gof_test(setosa, pnull, rnull, B=1000,
dnull=dnull, phat=phat, seed=seed0,
doMethods = "all")##> Using 23 cores.##>
##> MDgof multivariate goodness-of-fit tests
##> Data: continuous | n = 50 | dimension = 4 | B = 1012
##>
##> method statistic p.value
##> qKS 0.268 0.0000
##> qK 0.278 0.0010
##> qCvM 0.748 0.0000
##> qAD 8.419 0.1077
##> BR 1.099 0.7747
##>
##> Estimated parameters:
##> Sepal L. Sepal W. Petal L. Petal W.
##> 5.006000000 3.428000000 1.462000000 0.246000000 0.124248980
##>
##> 0.099216327 0.016355102 0.010330612 0.099216327 0.143689796
##>
##> 0.011697959 0.009297959 0.016355102 0.011697959 0.030159184
##>
##> 0.006069388 0.010330612 0.009297959 0.006069388 0.011106122MDgof::gof_test_adjusted_pvalue(setosa, pnull,
rnull, B=1000, dnull=dnull, phat=phat,
seed=seed0, doMethods = "all")##> Using 23 cores.##> p values of individual tests:##> qKS : 0##> qK : 0.001##> qCvM : 0##> qAD : 0.1077##> BR : 0.75##> adjusted p value of combined tests: 0##> qKS qK qCvM qAD BR
##> 0.0000000000 0.0009881423 0.0000000000 0.1077075099 0.7500000000
##> Min p
##> 0.0000000000The data can be found using
get_score_table <- function(url, max_goals = 4) {
goals <- read.csv(url)[, c("FTHG", "FTAG")]
with(
goals,
table(
factor(FTHG, levels = 0:max_goals),
factor(FTAG, levels = 0:max_goals)
)
)
}
urls <- c(
x = "https://www.football-data.co.uk/mmz4281/2021/E0.csv",
y = "https://www.football-data.co.uk/mmz4281/2526/E0.csv"
)
counts <- lapply(urls, get_score_table)
goals_premier_league <- as.matrix(cbind(
expand.grid(vals_x = 0:4, vals_y = 0:4),
lapply(counts, as.vector)))it is included here:
goals_premier_league <- cbind(
vals_x = c(0, 1, 2, 3, 4,
0, 1, 2, 3, 4,
0, 1, 2, 3, 4,
0, 1, 2, 3, 4,
0, 1, 2, 3, 4),
vals_y = c(rep(0:4, each = 5)),
x = c(30, 33, 26, 14, 6,
38, 40, 27, 13, 4,
23, 33, 7, 7, 3,
14, 18, 7, 6, 2,
3, 4, 3, 2, 0),
y = c(27, 22, 29, 27, 2,
31, 47, 37, 20, 6,
15, 27, 26, 11, 3,
10, 9, 10, 3, 2,
2, 4, 1, 2, 1)
)MD2sample::twosample_test(goals_premier_league,
seed=seed0, B=1000)##> Data is assumed to be discrete##> maxProcessor set to 1 because parallel overhead would dominate this short computation##>
##> MD2sample two-sample tests
##> Data: discrete | n.x = 363 | n.y = 374 | dimension = 2 | B = 1000
##>
##> method statistic p.value
##> KS 0.1314 0.0040
##> K 0.1827 0.0020
##> CvM 0.4884 0.0410
##> AD 2.3511 0.0530
##> NN 3.0350 0.0400
##> AZ 0.2111 0.0040
##> BF 0.3149 0.0040
##> ChiSquare 33.1600 0.0445MD2sample::twosample_test_adjusted_pvalue(goals_premier_league,
seed=seed0, B=1000, doMethods = "all")##> Data is assumed to be discrete
##> maxProcessor set to 1 because parallel overhead would dominate this short computation##> p values of individual tests:##> KS: 0.004##> K: 0.002##> CvM: 0.041##> AD: 0.053##> NN: 0.04##> AZ: 0.004##> BF: 0.004##> ChiSquare: 0.0445##> adjusted p value of combined tests: 0.013##> KS K CvM AD NN AZ
##> 0.00400000 0.00200000 0.04100000 0.05300000 0.04000000 0.00400000
##> BF ChiSquare Min p
##> 0.00400000 0.04446831 0.01300000