| 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 | # |