Commit | Line | Data |
---|---|---|
cbd88fe5 BA |
1 | context("OptimParams") |
2 | ||
3 | naive_f = function(link, M1,M2,M3, p,β,b) | |
4 | { | |
6dd5c2ac BA |
5 | d = length(M1) |
6 | K = length(p) | |
7 | λ <- sqrt(colSums(β^2)) | |
cbd88fe5 | 8 | |
6dd5c2ac BA |
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 | } | |
cbd88fe5 | 24 | |
6dd5c2ac BA |
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 | |
cbd88fe5 BA |
48 | } |
49 | ||
50 | test_that("naive computation provides the same result as vectorized computations", | |
51 | { | |
6dd5c2ac BA |
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] | |
cbd88fe5 | 58 | |
6dd5c2ac BA |
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)) | |
cbd88fe5 | 62 | |
6dd5c2ac BA |
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)) | |
cbd88fe5 | 67 | |
6dd5c2ac BA |
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) | |
cbd88fe5 | 75 | |
6dd5c2ac BA |
76 | # Test functions values |
77 | expect_equal( op$f(x), naive_f(link,M1,M2,M3, p,β,b) ) | |
cbd88fe5 | 78 | |
6dd5c2ac BA |
79 | # Test finite differences ~= gradient values |
80 | dir_h <- rep(0, (2+d)*K-1) | |
81 | dir_h[var] = h | |
cbd88fe5 | 82 | |
6dd5c2ac BA |
83 | expect_equal( op$grad_f(x)[var], ( op$f(x+dir_h) - op$f(x) ) / h, tol ) |
84 | } | |
85 | } | |
86 | } | |
cbd88fe5 | 87 | }) |
0f5fbd13 BA |
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 | # |