R/TestRankhist.R

Defines functions TestRankhist

Documented in TestRankhist

#' Statistical tests for rank histograms
#'
#' @description Perform statistical tests related to the deviation from flatness of a rank histogram.
#' @param rank.hist Vector of rank counts. Generated by function `Rankhist()`
#' @return A dataframe whose columns refer to the Pearson Chi^2 statistic, the Jolliffe-Primo test statistic for slope, and the Jolliffe-Primo test statistic for convexity. The rows refer to the actual test statistic and its p-value under the null hypothesis of a flat rank histogram.
#' @details
#' Given a vector of rank counts `x`, the Pearson Chi^2 statistic is calculated by
#' 
#' sum((x - sum(x)/length(x))^2 / (sum(x)/length(x)))
#' 
#' and has a chi^2 distribution with (length(x)-1) degrees of freedom if every rank is equally likely on average.
#' The Jolliffe-Primo test statistics are calculated by projecting the vector
#' 
#' (x-sum(x)/length(x)) / sqrt(sum(x)/length(x))
#' 
#' onto a linear, respectively squared contrast, i.e. a linear and quadratic function defined over the index set 1:length(x), who are mutually orthogonal, whose elements sum to zero, and whose squared elements sum to one. The projections independently have chi^2 distributions with 1 degree of freedom under the null hypothesis of a flat rank histogram.
#' @examples
#' data(eurotempforecast)
#' rh <- Rankhist(ens, obs)
#' TestRankhist(rh)
#' @seealso Rankhist, PlotRankhist
#' @references
#' Pearson K. (1900): X. On the criterion that a given system of deviations from the probable in the case of a correlated system of variables is such that it can be reasonably supposed to have arisen from random sampling.  Phil. Mag. Series 5, 50(302) \doi{10.1080/14786440009463897}
#' 
#' Jolliffe I.T., Primo C. (2008): Evaluating rank histograms using decompositions of the chi-square test statistic. Mon. Wea. Rev. 136(6) \doi{10.1175/2007MWR2219.1}
#' @export

TestRankhist <- function(rank.hist) {

  o.i <- rank.hist
  N <- sum(o.i)
  J <- length(o.i) 
  i <- 1:J
  e.i <- N / J
  x.i <- (o.i - e.i) / sqrt(e.i)
  # pearson chi^2 test 
  X2 <- sum(x.i * x.i)
  p.chisq <- pchisq(X2, df=J-1, lower.tail=FALSE)
  # jolliffe-primo: 
  # linear contrast
  a <- 2 * sqrt(3 / (J^3 - J))
  b <- -(sqrt(3) * J + sqrt(3)) / sqrt(J * (J + 1) * (J - 1))
  x.lin <- a*i+b
  X2.lin <- sum(x.i * x.lin)^2 # should have chi^2(df=1)
  # squared contrast
  a <- 6 * sqrt(5 / (J^5 - 5 * J^3 + 4 * J))
  b <- -1 / 2 * (sqrt(5) * J^2 - sqrt(5)) / 
       (sqrt((J - 2) * (J - 1) * J * (J + 1) * (J + 2)))
  x.u <- a * (i - (J + 1) / 2)^2 + b
  X2.u <- sum(x.i * x.u)^2 # should have chisq(df=1)

  # return 
  ret.df <- data.frame(pearson.chi2=c(X2, p.chisq), jp.lin=c(X2.lin, pchisq(X2.lin, df=1, lower.tail=FALSE)), jp.sq=c(X2.u, pchisq(X2.u, df=1, lower.tail=FALSE)))
  rownames(ret.df) <- c("test.statistic", "p.value")
  ret.df
}

Try the SpecsVerification package in your browser

Any scripts or data that you put into this service are public.

SpecsVerification documentation built on March 26, 2020, 7:55 p.m.