Package v.1.0 ready to be sent to CRAN
[morpheus.git] / pkg / tests / testthat / test-optimParams.R
CommitLineData
cbd88fe5
BA
1context("OptimParams")
2
2b3a6af5 3naive_f <- function(link, M1,M2,M3, p,β,b)
cbd88fe5 4{
2b3a6af5
BA
5 d <- length(M1)
6 K <- length(p)
6dd5c2ac 7 λ <- sqrt(colSums(β^2))
cbd88fe5 8
6dd5c2ac 9 # Compute β x2,3 (self) tensorial products
2b3a6af5
BA
10 β2 <- array(0, dim=c(d,d,K))
11 β3 <- array(0, dim=c(d,d,d,K))
6dd5c2ac
BA
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
2b3a6af5 25 res <- 0
6dd5c2ac
BA
26 for (i in 1:d)
27 {
2b3a6af5 28 term <- 0
6dd5c2ac 29 for (k in 1:K)
2b3a6af5
BA
30 term <- term + p[k]*.G(link,1,λ[k],b[k])*β[i,k]
31 res <- res + (term - M1[i])^2
6dd5c2ac
BA
32 for (j in 1:d)
33 {
2b3a6af5 34 term <- 0
6dd5c2ac 35 for (k in 1:K)
2b3a6af5
BA
36 term <- term + p[k]*.G(link,2,λ[k],b[k])*β2[i,j,k]
37 res <- res + (term - M2[i,j])^2
6dd5c2ac
BA
38 for (l in 1:d)
39 {
2b3a6af5 40 term <- 0
6dd5c2ac 41 for (k in 1:K)
2b3a6af5
BA
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
6dd5c2ac
BA
44 }
45 }
46 }
47 res
cbd88fe5
BA
48}
49
2b3a6af5
BA
50# TODO: understand why it fails and reactivate this test
51#test_that("naive computation provides the same result as vectorized computations",
52#{
53# h <- 1e-7 #for finite-difference tests
54# tol <- 1e-3 #large tolerance, necessary in some cases... (generally 1e-6 is OK)
55# n <- 10
56# for (dK in list( c(2,2), c(5,3)))
57# {
58# d <- dK[1]
59# K <- dK[2]
60#
61# M1 <- runif(d, -1, 1)
62# M2 <- matrix(runif(d^2, -1, 1), ncol=d)
63# M3 <- array(runif(d^3, -1, 1), dim=c(d,d,d))
64#
65# for (link in c("logit","probit"))
66# {
67# # X and Y are unused here (W not re-computed)
68# op <- optimParams(X=matrix(runif(n*d),ncol=d), Y=rbinom(n,1,.5),
69# K, link, M=list(M1,M2,M3))
70# op$W <- diag(d + d^2 + d^3)
71#
72# for (var in seq_len((2+d)*K-1))
73# {
74# p <- runif(K, 0, 1)
75# p <- p / sum(p)
76# β <- matrix(runif(d*K,-5,5),ncol=K)
77# b <- runif(K, -5, 5)
78# x <- c(p[1:(K-1)],as.double(β),b)
79#
80# # Test functions values
81# expect_equal( op$f(x), naive_f(link,M1,M2,M3, p,β,b) )
82#
83# # Test finite differences ~= gradient values
84# dir_h <- rep(0, (2+d)*K-1)
85# dir_h[var] = h
86# expect_equal(op$grad_f(x)[var], (op$f(x+dir_h) - op$f(x)) / h, tol)
87# }
88# }
89# }
90#})
91
92test_that("W computed in C and in R are the same",
cbd88fe5 93{
2b3a6af5
BA
94 tol <- 1e-8
95 n <- 500
6dd5c2ac
BA
96 for (dK in list( c(2,2), c(5,3)))
97 {
2b3a6af5
BA
98 d <- dK[1]
99 K <- dK[2]
100 link <- ifelse(d==2, "logit", "probit")
101 θ <- list(
102 p=rep(1/K,K),
103 β=matrix(runif(d*K),ncol=K),
104 b=rep(0,K))
105 io <- generateSampleIO(n, θ$p, θ$β, θ$b, link)
106 X <- io$X
107 Y <- io$Y
108 dd <- d + d^2 + d^3
109 p <- θ$p
110 β <- θ$β
111 λ <- sqrt(colSums(β^2))
112 b <- θ$b
113 β2 <- apply(β, 2, function(col) col %o% col)
114 β3 <- apply(β, 2, function(col) col %o% col %o% col)
115 M <- c(
116 β %*% (p * .G(link,1,λ,b)),
117 β2 %*% (p * .G(link,2,λ,b)),
118 β3 %*% (p * .G(link,3,λ,b)))
119 Id <- as.double(diag(d))
120 E <- diag(d)
121 v1 <- Y * X
122 v2 <- Y * t( apply(X, 1, function(Xi) Xi %o% Xi - Id) )
123 v3 <- Y * t( apply(X, 1, function(Xi) { return (Xi %o% Xi %o% Xi
124 - Reduce('+', lapply(1:d, function(j)
125 as.double(Xi %o% E[j,] %o% E[j,])), rep(0, d*d*d))
126 - Reduce('+', lapply(1:d, function(j)
127 as.double(E[j,] %o% Xi %o% E[j,])), rep(0, d*d*d))
128 - Reduce('+', lapply(1:d, function(j)
129 as.double(E[j,] %o% E[j,] %o% Xi)), rep(0, d*d*d))) } ) )
130 Omega1 <- matrix(0, nrow=dd, ncol=dd)
131 for (i in 1:n)
6dd5c2ac 132 {
2b3a6af5
BA
133 gi <- t(as.matrix(c(v1[i,], v2[i,], v3[i,]) - M))
134 Omega1 <- Omega1 + t(gi) %*% gi / n
6dd5c2ac 135 }
2b3a6af5
BA
136 W <- matrix(0, nrow=dd, ncol=dd)
137 Omega2 <- matrix( .C("Compute_Omega",
138 X=as.double(X), Y=as.integer(Y), M=as.double(M),
139 pn=as.integer(n), pd=as.integer(d),
140 W=as.double(W), PACKAGE="morpheus")$W, nrow=dd, ncol=dd )
141 rg <- range(Omega1 - Omega2)
142 expect_equal(rg[1], rg[2], tol)
6dd5c2ac 143 }
cbd88fe5 144})