pkg/tests/testthat/test-Moments_M3.R

   1 context("Moments_M3")
   2
   3 # (Slower, but trusted) R version of Moments_M3
   4 .Moments_M3_R = function(X,Y)
   5 {
   6   library(tensor)
   7   n = nrow(X)
   8   d = ncol(X)
   9   v = rep(0,d)
  10   e = diag(rep(1,d))
  11
  12   M31=M3_1=M3_2=M3_3 = tensor(tensor(v,v),v)
  13   for (i in 1:n)
  14     M31 = M31 + (Y[i]*tensor(tensor(X[i,],X[i,]),X[i,]))
  15   M31 = (1/n)*M31
  16
  17   for(j in 1:d)
  18   {
  19     l1=l2=l3 = tensor(tensor(v,v),v)
  20     for(i in 1:n)
  21     {
  22       l1 = l1 + Y[i]*tensor(tensor(e[,j],X[i,]),e[,j])
  23       l2 = l2 + Y[i]*tensor(tensor(e[,j],e[,j]),X[i,])
  24       l3 = l3 + Y[i]*tensor(tensor(X[i,],e[,j]),e[,j])
  25     }
  26     l1 = (1/n)*l1
  27     l2 = (1/n)*l2
  28     l3 = (1/n)*l3
  29     M3_1 = M3_1+l1
  30     M3_2 = M3_2+l2
  31     M3_3 = M3_3+l3
  32   }
  33
  34   M3 = M31-(M3_1+M3_2+M3_3)
  35   return (M3)
  36 }
  37
  38 test_that("both versions of Moments_M3 agree on various inputs",
  39 {
  40   for (n in c(20,200))
  41   {
  42     for (d in c(2,10,20))
  43     {
  44       X = matrix( runif(n*d,min=-1,max=1), nrow=n )
  45       Y = runif(n,min=-1,max=1)
  46       M3 = .Moments_M3(X,Y)
  47       M3_R = .Moments_M3_R(X,Y)
  48       expect_equal(max(abs(M3 - M3_R)), 0)
  49     }
  50   }
  51 })