R/boot.PC.r

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

################## PCBoot ##############################
#                                                      #
# Bootstraps PC model and returns an empirical         #
# empPC vector for the theta and r parameter.          #
########################################################

boot.PC <- function(model,size=100){

    if(class(model)!="PC") {
        stop("model has to be generated by PC");
    }
    
    P <-  attr(model,"P");
    cop <- attr(model,"cop");
    loss <- attr(model,"loss");
    domain <- attr(model,"domain");
    subdomains <- attr(model,"subdomains");
    
    sampleSize <- sum(P);
    Pnormalized <- P/sampleSize;

    thetaVector <- rep(NA,size);
    rVector <- rep(NA,size);
    lossVector <- rep(NA,size);
    sampleData <- list();

    for(i in 1:size){

        thisDraw <- rmultinom(1,sampleSize,Pnormalized);
        dim(thisDraw) <- dim(P);
        sampleData[[i]] <- thisDraw;

        thisResult <- PC(thisDraw,cop,loss,domain,subdomains=subdomains);
        thetaVector[i] <- attr(thisResult,"theta");
        rVector[i] <- attr(thisResult,"r");
        lossVector[i] <- attr(thisResult,"lossValue");
    }
    output <- list(model=model,theta=thetaVector,r=rVector,data=sampleData,loss=lossVector);
    class(output) <- "PCBoot";
    return(output);
}

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

    model <- x$model;
    rquantiles <- round(quantile(x$r,c(0.01,0.025,0.05,0.10,0.50,0.90,0.97,0.975,0.99)),3);
    tquantiles <- round(quantile(x$theta,c(0.01,0.025,0.05,0.10,0.50,0.90,0.97,0.975,0.99)),3);

    cat("\n--------- polyBoot Report ------------\n");
    cat("empPC Family:", attr(model,"cop"),"\n");
    cat("Loss Function Used:", attr(model,"loss"),"\n");
    cat("Estimated Theta:", attr(model,"theta"),"\n");
    cat("Estimated Correlation (r):", attr(model,"r"),"\n");
    cat("Standard Deviation for r:", sd(x$r),"\n");
    cat("\n------ Bootstrapped empPC ------\n");
    cat("Quantile\t Theta \t\t r\n");
    cat("0.01\t\t",rquantiles[1],"\t",tquantiles[1],"\n");
    cat("0.025\t\t",rquantiles[2],"\t",tquantiles[2],"\n");
    cat("0.05\t\t",rquantiles[3],"\t",tquantiles[3],"\n");
    cat("0.10\t\t",rquantiles[4],"\t",tquantiles[4],"\n");
    cat("0.50\t\t",rquantiles[5],"\t",tquantiles[5],"\n");
    cat("0.90\t\t",rquantiles[6],"\t",tquantiles[6],"\n");
    cat("0.95\t\t",rquantiles[7],"\t",tquantiles[7],"\n");
    cat("0.975\t\t",rquantiles[8],"\t",tquantiles[8],"\n");
    cat("0.99\t\t",rquantiles[9],"\t",tquantiles[9],"\n");
    cat("---------------------------------------\n");
}