tmp"
In MDgof: Various Methods for the Goodness-of-Fit Problem in D>1 Dimensions

knitr::opts_chunk$set(echo = TRUE, error=TRUE, 
                      warning = FALSE)
st=function() knitr::knit_exit()
library(mvtnorm)
library(MDgof)
library(copula)
B=100
M01=matrix(c(0, 1, 0, 1), 2, 2)
set.seed(111)

Binomial - No Estimation

pnull=function(x) {
  if(!is.matrix(x)) x=rbind(x)
  pbinom(x[,1], 5, 0.5)*pbinom(x[,2], 5, 0.6)
}
rnull=function() {
  nb=c(5,5)
  x=rbinom(1000, nb[1], 0.5)
  y=rbinom(1000, nb[2], 0.6)
  z=matrix(0,prod(nb+1),3)
  k=0
  for(i in 1:(nb[1]+1))
    for(j in 1:(nb[2]+1)) {
      k=k+1
      z[k, 1:2]=c(i,j)-1
      z[k, 3]=sum(x==i-1&y==j-1)
    }
  z
}  
ralt=function(p) {
  nb=c(5,5)
  x=rbinom(1000, nb[1], 0.5)
  y=rbinom(1000, nb[2], p)
  z=matrix(0,prod(nb+1),3)
  k=0
  for(i in 1:(nb[1]+1))
    for(j in 1:(nb[2]+1)) {
      k=k+1
      z[k, 1:2]=c(i,j)-1
      z[k, 3]=sum(x==i-1&y==j-1)
    }
  z
}

z=rnull()
gof_test(z, pnull, rnull,  B=B)
gof_power(pnull, rnull, ralt, c(0.6, 0.7), B=B)

Binomial - With Estimation

pnull=function(x, p) {
  if(!is.matrix(x)) x=rbind(x)
  pbinom(x[,1], 5, 0.5)*pbinom(x[,2], 5, p)
}
rnull=function(p) {
  nb=c(5,5)
  x=rbinom(1000, nb[1], 0.5)
  y=rbinom(1000, nb[2], p)
  z=matrix(0,prod(nb+1),3)
  k=0
  for(i in 1:(nb[1]+1))
    for(j in 1:(nb[2]+1)) {
      k=k+1
      z[k, 1:2]=c(i,j)-1
      z[k, 3]=sum(x==i-1&y==j-1)
    }
  z
}  
ralt=function(p) {
  nb=c(5,5)
  x=rbinom(1000, nb[1], p)
  y=rbinom(1000, nb[2], 0.6)
  z=matrix(0,prod(nb+1),3)
  k=0
  for(i in 1:(nb[1]+1))
    for(j in 1:(nb[2]+1)) {
      k=k+1
      z[k, 1:2]=c(i,j)-1
      z[k, 3]=sum(x==i-1&y==j-1)
    }
  z
}
phat=function(x) {
   nb=c(5, 5)
   z=aggregate(x[, 3], by = list(x[, 2]), FUN =sum)
   sum(z[,1]*z[,2])/nb[2]/sum(z[,2])

}

z=rnull(0.6)
gof_test(z, pnull, rnull,  phat=phat, B=B)
gof_power(pnull, rnull, rnull, c(0.5, 0.7), phat=phat, B=B)
gof_power(pnull, rnull, ralt, c(0.5, 0.7), phat=phat, B=B)

Correlated Normal - No Estimation

tmp=case.studies.disc(3, FALSE)
z=tmp$rnull()
gof_test(z, tmp$pnull, tmp$rnull,  B=B)
gof_power(tmp$pnull, tmp$rnull, tmp$ralt, 
          tmp$param_alt, B=B)

Correlated Normal - With Estimation

tmp=case.studies.disc(1, TRUE)
z=tmp$rnull(0)
gof_test(z, tmp$pnull, tmp$rnull,  phat=tmp$phat, B=B)
gof_power(tmp$pnull, tmp$rnull, tmp$ralt, 0.5, phat=tmp$phat, B=B)
z=tmp$ralt(5)
gof_test(z, tmp$pnull, tmp$rnull,  phat=tmp$phat, B=B)
gof_power(tmp$pnull, tmp$rnull, tmp$ralt, 8, phat=tmp$phat, B=B)