pkg/tests/testthat/test-Moments.R

   1 context("Moments")
   2
   3 test_that("both versions of Moments_Mi agree on various inputs",
   4 {
   5   for (n in c(20,200))
   6   {
   7     for (d in c(2,10,20))
   8     {
   9       E <- diag(d)
  10       Id <- as.double(E)
  11       X <- matrix( rnorm(n*d), nrow=n )
  12       Y <- rbinom(n, 1, .5)
  13       M2 <- as.double(.Moments_M2(X,Y))
  14       M2_R <- colMeans(Y * t( apply(X, 1, function(Xi) Xi %o% Xi - Id) ))
  15       expect_equal(max(abs(M2 - M2_R)), 0)
  16       M3 <- as.double(.Moments_M3(X,Y))
  17       M3_R <- colMeans(Y * t( apply(X, 1, function(Xi) {
  18         return (Xi %o% Xi %o% Xi -
  19           Reduce('+', lapply(1:d, function(j)
  20             as.double(Xi %o% E[j,] %o% E[j,])), rep(0, d*d*d)) -
  21           Reduce('+', lapply(1:d, function(j)
  22             as.double(E[j,] %o% Xi %o% E[j,])), rep(0, d*d*d)) -
  23           Reduce('+', lapply(1:d, function(j)
  24             as.double(E[j,] %o% E[j,] %o% Xi)), rep(0, d*d*d))) } ) ))
  25       expect_equal(max(abs(M3 - M3_R)), 0)
  26     }
  27   }
  28 })