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 <- 1e-3 #large tolerance, necessary in some cases... (generally 1e-6 is OK)
  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 })
  88
  89 # TODO: test computeW
  90 #    computeW = function(θ)
  91 #    {
  92 #      require(MASS)
  93 #      dd <- d + d^2 + d^3
  94 #      M <- Moments(θ)
  95 #      Id <- as.double(diag(d))
  96 #      E <- diag(d)
  97 #      v1 <- Y * X
  98 #      v2 <- Y * t( apply(X, 1, function(Xi) Xi %o% Xi - Id) )
  99 #      v3 <- Y * t( apply(X, 1, function(Xi) { return (Xi %o% Xi %o% Xi
 100 #        - Reduce('+', lapply(1:d, function(j) as.double(Xi %o% E[j,] %o% E[j,])), rep(0, d*d*d))
 101 #        - Reduce('+', lapply(1:d, function(j) as.double(E[j,] %o% Xi %o% E[j,])), rep(0, d*d*d))
 102 #        - Reduce('+', lapply(1:d, function(j) as.double(E[j,] %o% E[j,] %o% Xi)), rep(0, d*d*d))) } ) )
 103 #      Wtmp <- matrix(0, nrow=dd, ncol=dd)
 104 #
 105 #
 106 #g <- matrix(nrow=n, ncol=dd); for (i in 1:n) g[i,] = c(v1[i,], v2[i,], v3[i,]) - M
 107 #
 108 #
 109 #
 110 #
 111 #
 112 #
 113 #  p <- θ$p
 114 #  β <- θ$β
 115 #  b <- θ$b
 116 #
 117 #
 118 #
 119 #
 120 ##  # Random generation of the size of each population in X~Y (unordered)
 121 ##  classes <- rmultinom(1, n, p)
 122 ##
 123 ##  #d <- nrow(β)
 124 ##  zero_mean <- rep(0,d)
 125 ##  id_sigma <- diag(rep(1,d))
 126 ##  X <- matrix(nrow=0, ncol=d)
 127 ##  Y <- c()
 128 ##  for (i in 1:ncol(β)) #K = ncol(β)
 129 ##  {
 130 ##    newXblock <- MASS::mvrnorm(classes[i], zero_mean, id_sigma)
 131 ##    arg_link <- newXblock %*% β[,i] + b[i]
 132 ##    probas <-
 133 ##      if (li == "logit")
 134 ##      {
 135 ##        e_arg_link = exp(arg_link)
 136 ##        e_arg_link / (1 + e_arg_link)
 137 ##      }
 138 ##      else #"probit"
 139 ##        pnorm(arg_link)
 140 ##    probas[is.nan(probas)] <- 1 #overflow of exp(x)
 141 ##    X <- rbind(X, newXblock)
 142 ##    Y <- c( Y, vapply(probas, function(p) (rbinom(1,1,p)), 1) )
 143 ##  }
 144 #
 145 #
 146 #
 147 #
 148 #
 149 #
 150 #
 151 #
 152 #  Mhatt <- c(
 153 #    colMeans(Y * X),
 154 #    colMeans(Y * t( apply(X, 1, function(Xi) Xi %o% Xi - Id) )),
 155 #    colMeans(Y * t( apply(X, 1, function(Xi) { return (Xi %o% Xi %o% Xi
 156 #      - Reduce('+', lapply(1:d, function(j) as.double(Xi %o% E[j,] %o% E[j,])), rep(0, d*d*d))
 157 #      - Reduce('+', lapply(1:d, function(j) as.double(E[j,] %o% Xi %o% E[j,])), rep(0, d*d*d))
 158 #      - Reduce('+', lapply(1:d, function(j) as.double(E[j,] %o% E[j,] %o% Xi)), rep(0, d*d*d))) } ) ) ))
 159 #  λ <- sqrt(colSums(β^2))
 160 #  β2 <- apply(β, 2, function(col) col %o% col)
 161 #  β3 <- apply(β, 2, function(col) col %o% col %o% col)
 162 #  M <- c(
 163 #    β  %*% (p * .G(li,1,λ,b)),
 164 #    β2 %*% (p * .G(li,2,λ,b)),
 165 #    β3 %*% (p * .G(li,3,λ,b)) )
 166 #  print(sum(abs(Mhatt - M)))
 167 #
 168 #save(list=c("X", "Y"), file="v2.RData")
 169 #
 170 #
 171 #
 172 #
 173 #browser()
 174 #      for (i in 1:n)
 175 #      {
 176 #        gi <- t(as.matrix(c(v1[i,], v2[i,], v3[i,]) - M))
 177 #        Wtmp <- Wtmp + t(gi) %*% gi / n
 178 #      }
 179 #      Wtmp
 180 #      #MASS::ginv(Wtmp)
 181 #    },
 182 #
 183 #    #TODO: compare with R version?
 184 #    computeW_orig = function(θ)
 185 #    {
 186 #      require(MASS)
 187 #      dd <- d + d^2 + d^3
 188 #      M <- Moments(θ)
 189 #      Omega <- matrix( .C("Compute_Omega",
 190 #        X=as.double(X), Y=as.double(Y), M=as.double(M),
 191 #        pn=as.integer(n), pd=as.integer(d),
 192 #        W=as.double(W), PACKAGE="morpheus")$W, nrow=dd, ncol=dd )
 193 #      Omega
 194 #      #MASS::ginv(Omega) #, tol=1e-4)
 195 #    },
 196 #
 197 #    Moments = function(θ)
 198 #    {
 199 #      "Vector of moments, of size d+d^2+d^3"
 200 #
 201 #      p <- θ$p
 202 #      β <- θ$β
 203 #      λ <- sqrt(colSums(β^2))
 204 #      b <- θ$b
 205 #
 206 #      # Tensorial products β^2 = β2 and β^3 = β3 must be computed from current β1
 207 #      β2 <- apply(β, 2, function(col) col %o% col)
 208 #      β3 <- apply(β, 2, function(col) col %o% col %o% col)
 209 #
 210 #      c(
 211 #        β  %*% (p * .G(li,1,λ,b)),
 212 #        β2 %*% (p * .G(li,2,λ,b)),
 213 #        β3 %*% (p * .G(li,3,λ,b)))
 214 #    },
 215 #