pkg/tests/testthat/test-optimParams.R

   1 naive_f <- function(link, M1,M2,M3, p,β,b)
   2 {
   3   d <- length(M1)
   4   K <- length(p)
   5   λ <- sqrt(colSums(β^2))
   6
   7   # Compute β x2,3 (self) tensorial products
   8   β2 <- array(0, dim=c(d,d,K))
   9   β3 <- array(0, dim=c(d,d,d,K))
  10   for (k in 1:K)
  11   {
  12     for (i in 1:d)
  13     {
  14       for (j in 1:d)
  15       {
  16         β2[i,j,k] = β[i,k]*β[j,k]
  17         for (l in 1:d)
  18           β3[i,j,l,k] = β[i,k]*β[j,k]*β[l,k]
  19       }
  20     }
  21   }
  22
  23   res <- 0
  24   for (i in 1:d)
  25   {
  26     term <- 0
  27     for (k in 1:K)
  28       term <- term + p[k]*.G(link,1,λ[k],b[k])*β[i,k]
  29     res <- res + (term - M1[i])^2
  30     for (j in 1:d)
  31     {
  32       term <- 0
  33       for (k in 1:K)
  34         term <- term + p[k]*.G(link,2,λ[k],b[k])*β2[i,j,k]
  35       res <- res + (term - M2[i,j])^2
  36       for (l in 1:d)
  37       {
  38         term <- 0
  39         for (k in 1:K)
  40           term <- term + p[k]*.G(link,3,λ[k],b[k])*β3[i,j,l,k]
  41         res <- res + (term - M3[i,j,l])^2
  42       }
  43     }
  44   }
  45   res
  46 }
  47
  48 # TODO: understand why delta is so large (should be 10^-6 10^-7 ...)
  49 test_that("naive computation provides the same result as vectorized computations",
  50 {
  51   h <- 1e-7 #for finite-difference tests
  52   n <- 10
  53   for (dK in list( c(2,2), c(5,3)))
  54   {
  55     d <- dK[1]
  56     K <- dK[2]
  57
  58     M1 <- runif(d, -1, 1)
  59     M2 <- matrix(runif(d^2, -1, 1), ncol=d)
  60     M3 <- array(runif(d^3, -1, 1), dim=c(d,d,d))
  61
  62     for (link in c("logit","probit"))
  63     {
  64       # X and Y are unused here (W not re-computed)
  65       op <- optimParams(X=matrix(runif(n*d),ncol=d), Y=rbinom(n,1,.5),
  66         K, link, M=list(M1,M2,M3))
  67       op$W <- diag(d + d^2 + d^3)
  68
  69       for (var in seq_len((2+d)*K-1))
  70       {
  71         p <- runif(K, 0, 1)
  72         p <- p / sum(p)
  73         β <- matrix(runif(d*K,-5,5),ncol=K)
  74         b <- runif(K, -5, 5)
  75         x <- c(p[1:(K-1)],as.double(β),b)
  76
  77         # Test functions values (TODO: 1 is way too high)
  78         expect_equal( op$f(x)[1], naive_f(link,M1,M2,M3, p,β,b), tolerance=1 )
  79
  80         # Test finite differences ~= gradient values
  81         dir_h <- rep(0, (2+d)*K-1)
  82         dir_h[var] = h
  83         expect_equal( op$grad_f(x)[var], ((op$f(x+dir_h) - op$f(x)) / h)[1], tolerance=0.5 )
  84       }
  85     }
  86   }
  87 })
  88
  89 test_that("W computed in C and in R are the same",
  90 {
  91   tol <- 1e-8
  92   n <- 10
  93   for (dK in list( c(2,2))) #, c(5,3)))
  94   {
  95     d <- dK[1]
  96     K <- dK[2]
  97     link <- ifelse(d==2, "logit", "probit")
  98     θ <- list(
  99       p=rep(1/K,K),
 100       β=matrix(runif(d*K),ncol=K),
 101       b=rep(0,K))
 102     io <- generateSampleIO(n, θ$p, θ$β, θ$b, link)
 103     X <- io$X
 104     Y <- io$Y
 105     dd <- d + d^2 + d^3
 106     p <- θ$p
 107     β <- θ$β
 108     λ <- sqrt(colSums(β^2))
 109     b <- θ$b
 110     β2 <- apply(β, 2, function(col) col %o% col)
 111     β3 <- apply(β, 2, function(col) col %o% col %o% col)
 112     M <- c(
 113       β  %*% (p * .G(link,1,λ,b)),
 114       β2 %*% (p * .G(link,2,λ,b)),
 115       β3 %*% (p * .G(link,3,λ,b)))
 116     Id <- as.double(diag(d))
 117     E <- diag(d)
 118     v1 <- Y * X
 119     v2 <- Y * t( apply(X, 1, function(Xi) Xi %o% Xi - Id) )
 120     v3 <- Y * t( apply(X, 1, function(Xi) { return (Xi %o% Xi %o% Xi
 121       - Reduce('+', lapply(1:d, function(j)
 122         as.double(Xi %o% E[j,] %o% E[j,])), rep(0, d*d*d))
 123       - Reduce('+', lapply(1:d, function(j)
 124         as.double(E[j,] %o% Xi %o% E[j,])), rep(0, d*d*d))
 125       - Reduce('+', lapply(1:d, function(j)
 126         as.double(E[j,] %o% E[j,] %o% Xi)), rep(0, d*d*d))) } ) )
 127     Omega1 <- matrix(0, nrow=dd, ncol=dd)
 128     for (i in 1:n)
 129     {
 130       gi <- t(as.matrix(c(v1[i,], v2[i,], v3[i,]) - M))
 131       Omega1 <- Omega1 + t(gi) %*% gi / n
 132     }
 133     W <- matrix(0, nrow=dd, ncol=dd)
 134     Omega2 <- matrix( .C("Compute_Omega",
 135       X=as.double(X), Y=as.integer(Y), M=as.double(M),
 136       pnc=as.integer(1), pn=as.integer(n), pd=as.integer(d),
 137       W=as.double(W), PACKAGE="morpheus")$W, nrow=dd, ncol=dd )
 138     rg <- range(Omega1 - Omega2)
 139     expect_equal(rg[1], rg[2], tolerance=tol)
 140   }
 141 })