+test_that("W computed in C and in R are the same",
+{
+ tol <- 1e-8
+ n <- 500
+ for (dK in list( c(2,2), c(5,3)))
+ {
+ d <- dK[1]
+ K <- dK[2]
+ link <- ifelse(d==2, "logit", "probit")
+ θ <- list(
+ p=rep(1/K,K),
+ β=matrix(runif(d*K),ncol=K),
+ b=rep(0,K))
+ io <- generateSampleIO(n, θ$p, θ$β, θ$b, link)
+ X <- io$X
+ Y <- io$Y
+ dd <- d + d^2 + d^3
+ p <- θ$p
+ β <- θ$β
+ λ <- sqrt(colSums(β^2))
+ b <- θ$b
+ β2 <- apply(β, 2, function(col) col %o% col)
+ β3 <- apply(β, 2, function(col) col %o% col %o% col)
+ M <- c(
+ β %*% (p * .G(link,1,λ,b)),
+ β2 %*% (p * .G(link,2,λ,b)),
+ β3 %*% (p * .G(link,3,λ,b)))
+ Id <- as.double(diag(d))
+ E <- diag(d)
+ v1 <- Y * X
+ v2 <- Y * t( apply(X, 1, function(Xi) Xi %o% Xi - Id) )
+ v3 <- Y * t( apply(X, 1, function(Xi) { return (Xi %o% Xi %o% Xi
+ - Reduce('+', lapply(1:d, function(j)
+ as.double(Xi %o% E[j,] %o% E[j,])), rep(0, d*d*d))
+ - Reduce('+', lapply(1:d, function(j)
+ as.double(E[j,] %o% Xi %o% E[j,])), rep(0, d*d*d))
+ - Reduce('+', lapply(1:d, function(j)
+ as.double(E[j,] %o% E[j,] %o% Xi)), rep(0, d*d*d))) } ) )
+ Omega1 <- matrix(0, nrow=dd, ncol=dd)
+ for (i in 1:n)
+ {
+ gi <- t(as.matrix(c(v1[i,], v2[i,], v3[i,]) - M))
+ Omega1 <- Omega1 + t(gi) %*% gi / n
+ }
+ W <- matrix(0, nrow=dd, ncol=dd)
+ Omega2 <- matrix( .C("Compute_Omega",
+ X=as.double(X), Y=as.integer(Y), M=as.double(M),
+ pn=as.integer(n), pd=as.integer(d),
+ W=as.double(W), PACKAGE="morpheus")$W, nrow=dd, ncol=dd )
+ rg <- range(Omega1 - Omega2)
+ expect_equal(rg[1], rg[2], tol)
+ }