tests/testthat/test-score.R

In miheerdew/BMD-metrics: Metrics to compare cbce output.

context("test-score.R")
library(rlist)

set.seed(12345678)

n <- 70
dx <- 2
dy <- 3

epsilon <- 10 #SD of the noise term
A <- matrix(1, nrow=dx, ncol=dy)

test_that("Signal test for chisq score", {
  scores <- list.mapv(1:n,{
    X <- matrix(rnorm(dx*n), ncol=dx)
    Y <- X %*% A + matrix(rnorm(dy*n), ncol=dy)
    score(X, Y, method='chisq')
  })
  expect_lt(max(scores), 0.05)
})

test_that("Uniformity test for chisq score", {
  scores <- list.mapv(1:n, score(matrix(rnorm(dx*n), ncol=dx),
                                 matrix(rnorm(dy*n), ncol=dy),
                                 method='chisq'))
  expect_gt(ks.test(scores, y="punif")$p.value, .05)
})

test_that("Uniformity test for normal score", {
  scores <- list.mapv(1:n, score(matrix(rnorm(dx*n), ncol=dx),
                        matrix(rnorm(dy*n), ncol=dy),
                        method='normal'))
  expect_gt(ks.test(scores, y="punif")$p.value, .05)
})

test_that("Signal test for normal score", {
  scores <- list.mapv(1:n,{
              X <- matrix(rnorm(dx*n), ncol=dx)
              Y <- X %*% A + (sqrt(epsilon) * matrix(rnorm(dy*n), ncol=dy))
              score(X, Y, method='normal')
            })
  expect_lt(max(scores), 0.05)
})
test_that("Matrix is positive semidefinite", {
  S <- score(matrix(rnorm(dx*n), ncol=dx),
        matrix(rnorm(dy*n), ncol=dy), matrix_only=TRUE)$S
  list.mapv(1:n, {
    r <- rnorm(dx*dy)
    expect_gte(t(r) %*% S %*% r, 0)
  })
})

test_that("Matrices are correct", {
  X <- matrix(rnorm(dx*n), ncol=dx)
  Y <- matrix(rnorm(dy*n), ncol=dy)

  M <- score(X, Y, matrix_only=TRUE)
  N <- n - 1

  X <- scale(X)
  Y <- scale(Y)
  R <- cov(X, Y)

  from_pair <- function(x, y){
    #Encode x, y indices in single integer.
    x + dx*(y-1)
  }

  to_pair <- function(i) {
    #The inverse of from_pair
    x <- ((i-1) %% dx) + 1
    y <- (i - x)/dx + 1
    list(x=x, y=y)
  }

  # Check these indices for the matrices
  # Each index itself corresponds to a pair.

  indices <- list(c(1,1), c(1,2), c(2,1), c(2,2))

  list.iter(indices, ind ~ {
    p1 <- to_pair(ind[1])
    p2 <- to_pair(ind[2])

    expect_equal(M$S1[ind[1],ind[2]], sum(X[,p1$x]*X[,p2$x]*Y[,p1$y]*Y[,p2$y])/N)
    expect_equal(M$S2[ind[1], ind[2]],
      R[p1$x, p1$y]*R[p2$x,p2$y]*
        sum((X[,p1$x]^2 + Y[,p1$y]^2)*(X[,p2$x]^2 + Y[,p2$y]^2))/N
      )
    expect_equal(M$S3[ind[1], ind[2]],
            R[p1$x, p1$y]*sum((X[,p1$x]^2 + Y[,p1$y]^2)*X[,p2$x]*Y[,p2$y])/N)
  })
})