compare_chi_loss_funcs.R

## * Packages
library(mise)
mise()

## * Constnats .............................................
## Observed Values
obs <- rbind(
        c(1070, 1090, 1090),
        c(1090, 1020, 990)
    )
## Expected if H0 was true
e <- outer(rowSums(obs), colSums(obs)/sum(obs))
n = 10^5

## * Main Function. ............................................................

main <- function() {
    ## Built in Probability
    chiprob = chisq.test(obs)

    ## Print that probability
    print("Actual Probability:")
    print(chiprob$p.value)

    ## List all loss functions to test
    ## here loss function means difference from expected H0 distribution.
    loss_funcs <- list("RMSE" = loss_rmse,
                    "MAD" = loss_abs,
                    "Chi" = loss_chi)

    ## Calculate the p-value for different loss functions
    for (i in 1:length(loss_funcs)) {
        p <- pval_sim(obs, loss_funcs[[i]])
        print(names(loss_funcs)[i])
        print(p)
    }

    ## Compare them
    print(paste("Size:", sum(obs)))
    print(paste("Replicates", n))
}

## * Helper Functions...........................................................
expected_counts <- function(mat) {
    N <- sum(mat)
    e <- outer(rowSums(mat), colSums(mat)/N)
    return(e)
}

H0_sim <- function(mat) {
    multinomdata <- rmultinom(n = 1, prob = e, size = sum(obs))
    matrixData <- matrix(multinomdata, nrow = nrow(obs))
    return(matrixData)
}

## ** The P-value Replicate Simulation ........................................
pval_sim <- function(obs, loss) {
    ## Get the observed loss value (typicall chi=(o-e)^2/e)
    obs_loss <- loss(obs) 
    ## Get the Chi values if H0 was true
    chi_vals_null <- replicate(n, loss(H0_sim(obs)))
    ## Get the p-value by measuring the number of FalsePos
    my_chi_prob <- mean(chi_vals_null > obs_loss)

    return(my_chi_prob)
}

## ** Loss Functions ..........................................................
## ** Chi Squared .............................................................
loss_chi <- function(mat) {
    e <- expected_counts(mat)
    return(sum((mat-e)^2/e))
}
## ** RMSE ....................................................................
loss_rmse <- function(mat) {
    e     <- expected_counts(mat)
    se    <- (mat-e)^2
    mse   <- mean(se)  
    rmse  <- sqrt(mse)
    return(rmse)
}
## ** MAD .....................................................................
loss_abs <- function(mat) {
    e <- expected_counts(mat)
    return(sum(abs(mat-e)/e))
}

main()