gof.test <- function(X,E,fg="PowDiv",par=1,tpval="a",simrep=2000)

# funkce pro výpočet testu dobré shody

# vstupní parametry
#    X - vektor empirických četností
#    E - vektor teoretických četností
#    fg - typ testu, přípustné hodnoty jsou "PowDiv", "BWHD"
#     a "BWCS"
#    par - doplňující parametr, pro "PowDiv" může být libovolné
#     číslo, pro "BWHD" a "BWCS" číslo z intervalu (0,1)
#    tpval - způsob výpočtu p-hodnoty, přípustné hodnoty jsou
#     "a" (asymptotická p-hodnota počítaná z chi kvadrát rozdělení),
#     "e" (přesná p-hodnota, pouze pro malé vektory četností),
#     "s" (simulovaná p-hodnota počítaná metodou Monte Carlo)
#    simrep - počet opakování při výpočtu simulované p-hodnoty
# výstupní parametry
#    stat - testová statistika
#    df - počet stupňů volnosti
#    pval - p-hodnota testu

{
   name <- gof.name(fg,par);
   if(tpval=="e") {
      name <- paste(name,"with exact p-value",sep=" ");
   }
   if(tpval=="s") {
      name <- paste(name,"with simulated p-value (based on",
       simrep, "replicates)",sep=" ");
   }
   stat <- gof.stat(X,E,fg,par);
   df <- length(X)-1;
   pval <- gof.pval(X,E,fg,par,tpval,simrep);
   names(stat) <- "stat";
   names(df) <- "df";
   dname <- deparse(substitute(X));
   structure(list(statistic=stat,parameter=df,p.value=pval,
    method=name,data.name=dname,observed=X,expected=E),
    class="htest");
}


gof.name <- function(fg,par)

# vstupní parametry
#    fg - typ testu, přípustné hodnoty jsou "PowDiv", "BWHD"
#     a "BWCS"
#    par - doplňující parametr, pro "PowDiv" může být libovolné
#     číslo, pro "BWHD" a "BWCS" číslo z intervalu (0,1)
# výstupní parametry
#      methname - název statistiky

{
   if(((fg=="PowDiv")&&(par==1))||((fg=="BWHD")&&(par==0))
    ||((fg=="BWCS")&&(par==0))) {
      methname <- "Pearson chi-square test";
   }
   else if((fg=="PowDiv")&&(par==0)) {
      methname <- "Likelihood ratio test";
   }
   else if((fg=="PowDiv")&&(par==2/3)) {
      methname <- "Cressie-Read test";
   }
   else if((fg=="PowDiv")&&(par==-1/2)) {
      methname <- "Freeman-Tukey test";
   }
   else if((fg=="PowDiv")&&(par==-1)) {
      methname <- "Discriminant information test";
   }
   else if(((fg=="PowDiv")&&(par==-2))||((fg=="BWHD")&&(par==1))
    ||((fg=="BWCS")&&(par==1))) {
      methname <- "Neyman modified chi-square test";
   }
   else if(fg=="PowDiv") {
      methname <- paste("Power divergence test (lambda=",par,")",
       sep="");
   }
   else if(fg=="BWHD") {
      methname <- paste("Blended weight Hellinger distance test (",
       "alpha=",par,")",sep="")
   }
   else if(fg=="BWCS") {
      methname <- paste("Blended weight chi-square test (",
       "alpha=",par,")",sep="")
   }
   return(methname)
}


gof.stat <- function(X,E,fg,par,dig=NULL)

# funkce pro výpočet testové statistiky

# vstupní parametry
#    X - vektor empirických četností
#    E - vektor teoretických četností
#    fg - typ testu, přípustné hodnoty jsou "PowDiv", "BWHD"
#     a "BWCS"
#    par - doplňující parametr, pro "PowDiv" může být libovolné
#     číslo, pro "BWHD" a "BWCS" číslo z intervalu (0,1)
#    dig - počet desetinných míst při zaokrouhlování,
#     NULL pokud nemá zaokrouhlovat
# výstupní parametry
#      stat - testová statistika

{
   if(fg=="PowDiv") {
      if(par==0) { stat <- na.omit(2*X*log(X/E)) }
      else if(par==-1) { stat <-2*E*log(E/X) }
      else { stat <- (2*E*((X/E)^(par+1)-1))/(par*(par+1)) }
   }
   else if(fg=="BWHD") {
      stat <- ((X-E)^2)/((par*sqrt(X)+(1-par)*sqrt(E))^2)
   }
   else if(fg=="BWCS") {
      stat <- ((X-E)^2)/(par*X+(1-par)*E)
   }
   stat <- sum(stat)
   if(is.null(dig)==FALSE) { stat <- round(stat,dig) }
   return(stat)
}


gof.pval <- function(X=NULL,E,fg,par,tpval,simrep=2000,stat=NULL,
 stonly=FALSE,dig=NULL)

# funkce pro výpočet p-hodnoty testu

# vstupní parametry
#    X - vektor empirických četností
#    E - vektor teoretických četností
#    fg - typ testu, přípustné hodnoty jsou "PowDiv", "BWHD"
#     a "BWCS"
#    par - doplňující parametr, pro "PowDiv" může být libovolné
#     číslo, pro "BWHD" a "BWCS" číslo z intervalu (0,1)
#    tpval - způsob výpočtu p-hodnoty, přípustné hodnoty jsou
#     "a" (asymptotická p-hodnota počítaná z chi kvadrát rozdělení),
#     "e" (přesná p-hodnota, pouze pro malé vektory četností),
#     "s" (simulovaná p-hodnota počítaná metodou Monte Carlo)
#    simrep - počet opakování při výpočtu simulované p-hodnoty
#    stat - statistika, volitelné, jinak se vypočítá z X
#    stonly - TRUE pokud má místo p-hodnoty vypočítat
#      pouze pravděpodobnost statistiky
#    dig - počet desetinných míst při zaokrouhlování,
#     NULL pokud nemá zaokrouhlovat
# výstupní parametry
#      pvalue - p-hodnota testu

#   používá funkci rmultz2 z balíku combinat

{
   digst <- 10;
   k <- length(E);
   n <- sum(E);
   P <- E/n;
   if(is.null(stat)) { stat <- gof.stat(X,E,fg,par,dig=digst); }
   if(tpval=="a") {
      pvalue <- 1-pchisq(stat,k-1);
   }
   else if(tpval=="e") {
      mat <- matAll(n,k);
      if(stonly==FALSE) {
         mat <- matrix(mat[apply(mat,MARGIN=1,FUN=gof.stat,E=E,
          fg=fg,par=par,dig=digst)>=stat,],ncol=k);
      } else {
         mat <- matrix(mat[apply(mat,MARGIN=1,FUN=gof.stat,E=E,
          fg=fg,par=par,dig=digst)==stat,],ncol=k);
      }
      pvalue <- sum(apply(mat,MARGIN=1,FUN=dmultinom,prob=P))
   }
   else if(tpval=="s") {
      mat <- t(rmultz2(n,P,simrep));
      matrow <- nrow(mat);
      if(stonly==FALSE) {
         mat <- matrix(mat[apply(mat,MARGIN=1,FUN=gof.stat,E=E,
          fg=fg,par=par,dig=digst)>=stat,],ncol=k);
      } else {
         mat <- matrix(mat[apply(mat,MARGIN=1,FUN=gof.stat,E=E,
          fg=fg,par=par,dig=digst)==stat,],ncol=k);
      }
      pvalue <- nrow(mat)/matrow;
   }
   if(is.null(dig)==FALSE) { pvalue <- round(pvalue,dig) }
   return(pvalue);
}


gof.power <- function(n,P0,P1,size,fg,par,tpval="e",rand=FALSE,
 dig=NULL)

# funkce pro výpočet pravděpodobosti zamítnutí hypotézy P0
#  za platnosti P1, pouze pro malé vektory četností

# vstupní parametry
#    n - součet hodnot vektoru četností
#    P0 - vektor pravděpodobností nulové hypotézy
#    P1 - vektor pravděpodobností alternativní hypotézy
#    size - přesná hladina testu
#    fg - typ testu, přípustné hodnoty jsou "PowDiv", "BWHD"
#     a "BWCS"
#    par - doplňující parametr, pro "PowDiv" může být libovolné
#     číslo, pro "BWHD" a "BWCS" číslo z intervalu (0,1)
#    tpval - způsob výpočtu p-hodnoty, přípustné hodnoty jsou
#     "a" (asymptotická p-hodnota počítaná z chi kvadrát rozdělení),
#     "e" (přesná p-hodnota, pouze pro malé vektory četností),
#     "s" (simulovaná p-hodnota počítaná metodou Monte Carlo)
#    rand - TRUE pokud má použít znáhodněný test
#    dig - počet desetinných míst při zaokrouhlování,
#     NULL pokud nemá zaokrouhlovat
# výstupní parametry
#    power - síla testu

{
   digst <- 10;
   k <- length(P0);
   E0 <- P0*n;
   mat <- matAll(n,k);
   maxStat <- max(apply(matrix(mat[apply(mat,MARGIN=1,
    FUN=gof.pval,E=E0,fg=fg,par=par,tpval=tpval)>=size,],
    ncol=k),MARGIN=1,FUN=gof.stat,E=E0,fg=fg,par=par,
    dig=digst));
   mat1 <- matrix(mat[apply(mat,MARGIN=1,FUN=gof.stat,E=E0,
    fg=fg,par=par,dig=digst)>maxStat,],ncol=k);
   mat2 <- matrix(mat[apply(mat,MARGIN=1,FUN=gof.stat,E=E0,
    fg=fg,par=par,dig=digst)==maxStat,],ncol=k);
   if(rand) {
      pran <- (size-gof.pval(E=E0,fg=fg,par=par,tpval=tpval,
       stat=maxStat))/gof.pval(E=E0,fg=fg,par=par,tpval=tpval,
       stat=maxStat,stonly=TRUE)+1;
   } else {
      pran <- 0;
   }
   power <- (sum(apply(mat1,MARGIN=1,FUN=dmultinom,prob=P1))
    +pran*sum(apply(mat2,MARGIN=1,FUN=dmultinom,prob=P1)));
   if(is.null(dig)==FALSE) { power <- round(power,dig); }
   return(power);
}


matAll <- function(n,k)

# funkce pro výpočet matice všech možných kombinací
#  vektoru četností

# vstupní parametry
#    n - součet hodnot vektoru četností
#    k - velikost vektoru četností
# výstupní parametry
#    mat - matice možností

{
   egA <- 0:n;
   egB <- list(egA);
   egC <- egB;
   for(i in 1:(k-1)) { egC <- c(egC,egB); }
   mat <- as.matrix(expand.grid(egC));
   mat <- mat[rowSums(mat)==n,];
   mat <- matrix(mat,ncol=k);
   return(mat);
}