R/GF.r

In jeroenooms/JJcorr: Calculates polychorical correlations for several copula families.

################ GF Goodness of Fit #################
#                                                         #
# This function performs a Goodness of Fit test using     #
# simulation. It resamples data from the fitted model and #
# calculates the loss of the simulated data with respect  #
# to the model. The loss-empPC is then compared    #
# with the loss of the observed data with respect to the  #
# model.                                                  #
#                                                         #
# By default, the same loss function is used to evaluate  #
# GF as is used to fit the model. However this is not     #
# required.                                               #
###########################################################

GF <- function(model,size=10000,fitLoss="modelLoss"){

    if(fitLoss=="modelLoss") fitLoss <- attr(model,"loss");

    if(class(model)!="PC") {
        stop("model has to be generated by PC");
    }

    P <-  attr(model,"P");
    loss <- attr(model,"loss");

    fittedTable <- attr(model,"fitted");
    sampleSize <- sum(P);

    pNormalized <- P/sum(P);
    fNormalized <- fittedTable/sum(fittedTable);
    fitLossFunction <- get(paste("loss",fitLoss,sep=""));
    
    diffVector <- rep(NA,size);

    for(i in 1:size){

        thisDraw <- rmultinom(1,sampleSize,fNormalized);
        dim(thisDraw) <- dim(P);
        drawNormalized <- thisDraw/sum(thisDraw);
        
        diffVector[i] <- fitLossFunction(fittedTable,drawNormalized);
    }
    modelLoss <- fitLossFunction(fittedTable,pNormalized);
    p <- mean(modelLoss < diffVector)
    
    output <- list(model=model,fitLossValue=modelLoss,fitLossDist=diffVector,p=p);
    attr(output,"fitLoss") <- fitLoss;
    class(output) <- "GF";
    return(output);
}

print.GF <- function(x, ...){

    print(x$model);
    
    cat("Goodness of Fit:\n");
    cat("fit loss function used:",attr(x,"fitLoss"),"\n");
    cat("Observed loss:",x$fitLossValue,", p=",x$p,"\n");
    cat("---------------------------------------\n");
}