Simplify data sampling + cosmetics + draft unit test for W matrix in optim
[morpheus.git] / 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 #