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)
})
})
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.