poissonGoodnessOfFit <- function(BSseq, nQuantiles = 10^5) {
    Cov <- getBSseq(BSseq, type = "Cov")
    lambda <- rowMeans2(Cov)
    out <- list(chisq = NULL, df = ncol(Cov) - 1)
    out$chisq <- rowSums2((Cov - lambda)^2 / lambda)
    if(nQuantiles > length(out$chisq))
        nQuantiles <- length(out$chisq)
    probs <- seq(from = 0, to = 1, length.out = nQuantiles)
    out$quantiles <- quantile(out$chisq, prob = probs, na.rm = TRUE)
    class(out) <- "chisqGoodnessOfFit"
    return(out)
}

print.chisqGoodnessOfFit <- function(x, ...) {
    cat("Chisq Goodness of Fit object\n")
    cat(" ", length(x$chisq), "tests\n")
    cat(" ", x$df, "degrees of freedom\n")
}

binomialGoodnessOfFit <- function(BSseq, method = c("MLE"), nQuantiles = 10^5) {
    method <- match.arg(method)
    M <- getCoverage(BSseq, type = "M", what = "perBase")
    Cov <- getCoverage(BSseq, type = "Cov", what = "perBase")
    if(method == "MLE") # This is the MLE under iid
        p <- rowSums2(M) / rowSums2(Cov)
    ## p <- rowMeans(M / Cov, na.rm = TRUE)
    out <- list(chisq = NULL, quantiles = NULL, df = ncol(Cov) - 1)
    ## This gets tricky.  For the binomial model, we need Cov > 0, but
    ## this implies that a given location may have less then nCol(BSseq)
    ## number of observations, and hence the degrees of freedom will
    ## vary.
    out$chisq <- rowSums2((M - Cov*p)^2 / sqrt(Cov * p * (1-p)))
    probs <- seq(from = 0, to = 1, length.out = nQuantiles)
    ## Next line will remove all locations where not all samples have coverage
    out$quantiles <- quantile(out$chisq, prob = probs, na.rm = TRUE)
    class(out) <- "chisqGoodnessOfFit"
    return(out)
}

plot.chisqGoodnessOfFit <- function(x, type = c("chisq", "pvalue"),
                                    plotCol = TRUE, qqline = TRUE, pch = 16, cex = 0.75, ...) {
    type <- match.arg(type)
    switch(type,
           "chisq" = {
               yy <- x$quantiles
               xx <- qchisq(ppoints(yy), df = x$df)
           },
           "pvalue" = {
               yy <- sort(1 - qchisq(x$quantiles, df = x$df))
               xx <- qunif(ppoints(yy), min = 0, max = 1)
           })
    if(plotCol) {
        colors <- rep("black", length(yy))
        colors[floor(length(colors)*.95):length(colors)] <- "red"
        colors[floor(length(colors)*.99):length(colors)] <- "violet"
        colors[floor(length(colors)*.999):length(colors)] <- "orange"
    } else {
        colors <- "black"
    }
    qqplot(xx, yy, xlab = "Theoretical quantiles",
           ylab = "Observed quantiles", col = colors, pch = pch, cex = cex, ...)
    if(qqline) {
        yyqt <- quantile(yy, c(0.25, 0.75))
        xxqt <- quantile(xx, c(0.25, 0.75))
        slope <- diff(yyqt) / diff(xxqt)
        int <- yyqt[1L] - slope * xxqt[1L]
        abline(int, slope, col = "blue")
    }
    abline(0, 1, col = "blue", lty = 2)
}