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