| 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 | }) |