pkg/tests/testthat/test-optimParams.R

   1 context("OptimParams")
   2
   3 naive_f = function(link, M1,M2,M3, p,β,b)
   4 {
   5     d = length(M1)
   6     K = length(p)
   7     λ <- sqrt(colSums(β^2))
   8
   9     # Compute β x2,3 (self) tensorial products
  10     β2 = array(0, dim=c(d,d,K))
  11     β3 = array(0, dim=c(d,d,d,K))
  12     for (k in 1:K)
  13     {
  14         for (i in 1:d)
  15         {
  16             for (j in 1:d)
  17             {
  18                 β2[i,j,k] = β[i,k]*β[j,k]
  19                 for (l in 1:d)
  20                     β3[i,j,l,k] = β[i,k]*β[j,k]*β[l,k]
  21             }
  22         }
  23     }
  24
  25     res = 0
  26     for (i in 1:d)
  27     {
  28         term = 0
  29         for (k in 1:K)
  30             term = term + p[k]*.G(link,1,λ[k],b[k])*β[i,k]
  31         res = res + (term - M1[i])^2
  32         for (j in 1:d)
  33         {
  34             term = 0
  35             for (k in 1:K)
  36                 term = term + p[k]*.G(link,2,λ[k],b[k])*β2[i,j,k]
  37             res = res + (term - M2[i,j])^2
  38             for (l in 1:d)
  39             {
  40                 term = 0
  41                 for (k in 1:K)
  42                     term = term + p[k]*.G(link,3,λ[k],b[k])*β3[i,j,l,k]
  43                 res = res + (term - M3[i,j,l])^2
  44             }
  45         }
  46     }
  47     res
  48 }
  49
  50 test_that("naive computation provides the same result as vectorized computations",
  51 {
  52     h <- 1e-7 #for finite-difference tests
  53     tol <- 5e-4 #.25 * sqrt(h) #about 7.9 e-5
  54     for (dK in list( c(2,2), c(5,3)))
  55     {
  56         d = dK[1]
  57         K = dK[2]
  58
  59         M1 = runif(d, -1, 1)
  60         M2 = matrix(runif(d*d,-1,1), ncol=d)
  61         M3 = array(runif(d*d*d,-1,1), dim=c(d,d,d))
  62
  63         for (link in c("logit","probit"))
  64         {
  65             op = new("OptimParams", "li"=link, "M1"=as.double(M1),
  66                 "M2"=as.double(M2), "M3"=as.double(M3), "K"=as.integer(K))
  67
  68             for (var in seq_len((2+d)*K-1))
  69             {
  70                 p = runif(K, 0, 1)
  71                 p = p / sum(p)
  72                 β <- matrix(runif(d*K,-5,5),ncol=K)
  73                 b = runif(K, -5, 5)
  74                 x <- c(p[1:(K-1)],as.double(β),b)
  75
  76                 # Test functions values
  77                 expect_equal( op$f(x), naive_f(link,M1,M2,M3, p,β,b) )
  78
  79                 # Test finite differences ~= gradient values
  80                 dir_h <- rep(0, (2+d)*K-1)
  81                 dir_h[var] = h
  82
  83                 expect_equal( op$grad_f(x)[var], ( op$f(x+dir_h) - op$f(x) ) / h, tol )
  84             }
  85         }
  86     }
  87 })