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

knitr::opts_chunk$set(echo = TRUE, error=TRUE)
library(mvtnorm)
library(MDgof)
st=function() knitr::knit_exit()
set.seed(113)
M01=matrix(c(0, 1, 0, 1), 2, 2)
B=50

Continuous Data - Dim=2 - No Estimation

pnullU=function(x) {
  if(!is.matrix(x)) x=rbind(x)
  d=ncol(x)
  y=rep(0, nrow(x))
  for(i in seq_along(y)) {
    y[i]=punif(x[i,1])*punif(x[i,2])
    if(d>2) 
      for(j in 3:d) y[i]=y[i]*punif(x[i,j])
  }  
  y
}
dnullU=function(x) {
  if(!is.matrix(x)) return(1) 
  rep(1, nrow(x))
}  
rnullU=function() round(matrix(runif(2*250), ncol=2), 3)
raltU=function(a=1) round(matrix(rbeta(2*250, a, a), ncol=2), 3)

pnullN=function(x) {
  if(!is.matrix(x)) x=rbind(x)
  y=rep(0, nrow(x))
  for(i in seq_along(y))
    y[i]=mvtnorm::pmvnorm(rep(-Inf, ncol(x)), x[i, ])
  y
}
dnullN=function(x) {
  if(!is.matrix(x)) x=rbind(x)
  k=nrow(x)
  y=rep(0, k)
  for(i in 1:k)
      y[i]=mvtnorm::dmvnorm(x[i,], mean=rep(0,ncol(x))) 
  y
}
rnullN=function() round(mvtnorm::rmvnorm(250, c(0,0)), 3)
raltN=function(a=0) round(mvtnorm::rmvnorm(250, c(a, a)), 3)

dta=rnullU()
gof_test(dta, pnullU, rnullU, Ranges=M01, dnull=dnullU, B=B)
dta=raltU(1.5)
gof_test(dta, pnullU, rnullU, Ranges=M01, dnull=dnullU, B=B)
gof_test(dta, pnullU, rnullU, TS=TS_cont, dnull=dnullU, B=B)
dta=rnullN()
gof_test(dta, pnullN, rnullN, dnull=dnullU, B=B)
gof_power(pnullU, rnullU, raltU, c(1, 1.4), Ranges=M01, 
          dnull=dnullU, B=B)
gof_power(pnullN, rnullN, raltN, c(0, 0.25), 
          dnull=dnullN, B=B)
gof_power(pnullU, rnullU, raltU, c(1, 1.5), 
          dnull=dnullU, TS=TS_cont,   B=B)

Continuous Data - Dim=5 - No Estimation

pnullU=function(x) {
  if(!is.matrix(x)) x=rbind(x)
  d=ncol(x)
  y=rep(0, nrow(x))
  for(i in seq_along(y)) {
    y[i]=punif(x[i,1])*punif(x[i,2])
    if(d>2) 
      for(j in 3:d) y[i]=y[i]*punif(x[i,j])
  }  
  y
}
dnullU=function(x) {
  if(!is.matrix(x)) return(1)
  rep(1, nrow(x))
}
raltU=function(a=1) round(matrix(rbeta(5*250, a, a), ncol=5), 3)
rnullU=function() round(matrix(runif(5*250), ncol=5), 3)

gof_test(rnullU(), pnullU, rnullU, dnull=dnullU, B=B)
gof_test(rnullU(), pnullU, rnullU, TS=TS_cont, 
         dnull=dnullU, B=B)
gof_power(pnullU, rnullU, raltU, c(1, 1.5),  
          dnull=dnullU, B=B)
gof_power(pnullU, rnullU, raltU, c(1, 1.5),  
          dnull=dnullU, TS=TS_cont, B=B)

Continuous - Dim=2 - with Estimation

pnullNE=function(x, mu=rep(0,ncol(x))) {
  if(!is.matrix(x)) x=rbind(x)
  k=nrow(x)
  y=rep(0, k)
  for(i in 1:k)
      y[i]=mvtnorm::pmvnorm(lower = rep(-Inf, ncol(x)),
               upper = x[i,], mean  = mu, sigma=diag(ncol(x))) 
  y
}  
rnullNE=function(mu=c(0,0)) mvtnorm::rmvnorm(250, mean=mu)
ralt_mu=function(mu) mvtnorm::rmvnorm(250, mean=c(mu,mu))
ralt_sig=function(s=0) 
    mvtnorm::rmvnorm(250, sigma = matrix(c(1,s,s,1),2,2))
dnullNE=function(x, mu=rep(0,ncol(x))) {
  if(!is.matrix(x)) x=rbind(x)
  k=nrow(x)
  y=rep(0, k)
  for(i in 1:k)
      y[i]=mvtnorm::dmvnorm(x[i,], mean  = mu) 
  y
}
phat=function(x) apply(x,2,mean)

gof_test(ralt_mu(0), pnullNE, rnullNE, phat=phat, 
         dnull=dnullNE, B=B)
gof_test(ralt_mu(1), pnullNE, rnullNE, phat=phat, 
         dnull=dnullNE, TS=TS_cont, B=B)
gof_test(ralt_sig(0.6), pnullNE, rnullNE, phat=phat, 
         dnull=dnullNE, B=B)
gof_power(pnullNE, rnullNE, ralt_mu, c(0, 0.5), 
          dnull=dnullNE, phat=phat, B=B)
gof_power(pnullNE, rnullNE, ralt_sig, c(0, 0.5), 
          dnull=dnullNE, phat=phat, B=B)

pnull=function(x, s) {
         f=function(x) mvtnorm::pmvnorm(rep(-Inf, length(x)), x,
                        mean=s$mu,                
                        sigma=s$sigma)
       if(!is.matrix(x)) return(f(x))
       apply(x, 1, f)      
}
dnull=function(x, s) {
         f=function(x) mvtnorm::dmvnorm(x,
                        mean=s$mu,                
                        sigma=s$sigma)
       if(!is.matrix(x)) return(f(x))
       apply(x, 1, f)      
}
rnull = function(s=list(mu=rep(0, 2),sigma=diag(2))) {
  mvtnorm::rmvnorm(100, mean=s$mu, sigma=s$sigma)
}
ralt=function(s=0) 
    mvtnorm::rmvnorm(100, sigma = matrix(c(1,s,s,1),2,2))
phat=function(x) {list(mu=apply(x, 2, mean), sigma=cov(x))}      
x=rnull()
MDgof::gof_test(x, pnull, rnull, phat=phat, B=50)
gof_power(pnull, rnull, ralt, c(0, 0.5), 
          dnull=dnull, phat=phat, B=B)