set.seed(290875)

library("party")

### get rid of the NAMESPACE

attach(asNamespace("party"))

### 

###

###    Regression tests for linear statistics, expectations and covariances

###    

###    functions defined in file `./src/LinearStatistics.c'    

### tests for function C_LinearStatistic

### Linear Statistics

x = matrix(c(rep.int(1,4), rep.int(0,6)), ncol = 1)

y = matrix(1:10, ncol = 1)

weights = rep(1, 10)

linstat = LinearStatistic(x, y, weights)

stopifnot(isequal(linstat, sum(1:4)))

weights[1] = 0

linstat = LinearStatistic(x, y, weights) 

stopifnot(isequal(linstat, sum(2:4)))

xf <- gl(3, 10)

yf <- gl(3, 10)[sample(1:30)]

x <- sapply(levels(xf), function(l) as.numeric(xf == l))

colnames(x) <- NULL

y <- sapply(levels(yf), function(l) as.numeric(yf == l))

colnames(y) <- NULL

weights <- sample(1:30)

linstat <- LinearStatistic(x, y, weights) 

stopifnot(isequal(linstat, as.vector(t(x) %*% diag(weights) %*% y)))

xf <- factor(cut(rnorm(6000), breaks = c(-Inf, -2, 0.5, Inf)))

x <- sapply(levels(xf), function(l) as.numeric(xf == l))

yf <- factor(cut(rnorm(6000), breaks = c(-Inf, -0.5, 1.5, Inf)))

y <- sapply(levels(yf), function(l) as.numeric(yf == l))

weights <- rep(1, nrow(x))

colnames(x) <- NULL

colnames(y) <- NULL

weights <- rep(1, 6000)

linstat <- LinearStatistic(x, y, weights)

stopifnot(isequal(as.vector(table(xf, yf)), linstat))

stopifnot(isequal(as.vector(t(x)%*%y), linstat))

### tests for function C_ExpectCovarInfluence

eci <- ExpectCovarInfluence(y, weights)

isequal(eci@sumweights, sum(weights))

isequal(eci@expectation, drop(weights %*% y / sum(weights)))

ys <- t(t(y) - eci@expectation)

stopifnot(isequal(eci@covariance, (t(ys) %*% (weights * ys)) /

                  sum(weights)))

### tests for function C_ExpectCovarLinearStatistic

### Conditional Expectation and Variance (via Kruskal-Wallis statistic)

### case 1: p > 1, q = 1

group <- gl(3, 5)

x <- sapply(levels(group), function(l) as.numeric(group == l))

y <- matrix(1:15, ncol = 1)

weights <- rep(1, 15)

linstat <- LinearStatistic(x, y, weights)

expcov <- ExpectCovarLinearStatistic(x, y, weights)

KW <- quadformTestStatistic(linstat, expcov@expectation, expcov@covariance)

kts <- kruskal.test(y ~ group)$statistic

stopifnot(isequal(KW, kts))

### case 2: p = 1, q > 1

linstat <- LinearStatistic(y, x, weights)

expcov <- ExpectCovarLinearStatistic(y, x, weights)

KW <- quadformTestStatistic(linstat, expcov@expectation, expcov@covariance)

kts <- kruskal.test(y ~ group)$statistic

stopifnot(isequal(KW, kts))

### case 3: p = 1, q = 1

x <- x[,1,drop = FALSE]

linstat <- LinearStatistic(x, y, weights)

expcov <- ExpectCovarLinearStatistic(x, y, weights)

KW <- quadformTestStatistic(linstat, expcov@expectation, expcov@covariance)

kts <- kruskal.test(y ~ as.factor(x))$statistic

stopifnot(isequal(KW, kts))

### case 4: p > 1, q > 1 via chisq.test

n <- 900

xf <- gl(3, n / 3)

yf <- gl(3, n / 3)[sample(1:n)]

x <- sapply(levels(xf), function(l) as.numeric(xf == l))

colnames(x) <- NULL

y <- sapply(levels(yf), function(l) as.numeric(yf == l))

colnames(y) <- NULL

weights <- rep(1, n)

linstat <- LinearStatistic(x, y, weights)

expcov <- ExpectCovarLinearStatistic(x, y, weights)

chi <- quadformTestStatistic(linstat, expcov@expectation, expcov@covariance)

chis <- chisq.test(table(xf, yf))$statistic

stopifnot(isequal(round(chi, 1), round(chis, 1)))

### tests for function C_PermutedLinearStatistic

### Linear Statistics with permuted indices

x <- matrix(rnorm(100), ncol = 2)

y <- matrix(rnorm(100), ncol = 2)

weights <- rep(1, 50)

indx <- 1:50

perm <- 1:50

stopifnot(isequal(LinearStatistic(x, y, weights), 

                  PermutedLinearStatistic(x, y, indx, perm)))

x <- matrix(1:10000, ncol = 2)

y <- matrix(1:10000, ncol = 2)

for (i in 1:100) {

  indx <- sample(1:ncol(y), replace = TRUE)

  perm <- sample(1:ncol(y), replace = TRUE)

  stopifnot(isequal(as.vector(t(x[indx,]) %*% y[perm, ]),

                    PermutedLinearStatistic(x, y, indx, perm)))

}