#' o$f( o$linArgs(par0) )
#' o$f( o$linArgs(par1) )
#' @export
-optimParams = function(X, Y, K, link=c("logit","probit"))
+optimParams <- function(X, Y, K, link=c("logit","probit"))
{
# Check arguments
if (!is.matrix(X) || any(is.na(X)))
stop("X: numeric matrix, no NAs")
- if (!is.numeric(Y) || any(is.na(Y)) || any(Y!=0 | Y!=1))
+ if (!is.numeric(Y) || any(is.na(Y)) || any(Y!=0 & Y!=1))
stop("Y: binary vector with 0 and 1 only")
link <- match.arg(link)
if (!is.numeric(K) || K!=floor(K) || K < 2)
stop("Missing arguments")
# Precompute empirical moments
- M <- computeMoments(optargs$X,optargs$Y)
+ M <- computeMoments(X, Y)
M1 <- as.double(M[[1]])
M2 <- as.double(M[[2]])
M3 <- as.double(M[[3]])
- Mhat <<- matrix(c(M1,M2,M3), ncol=1)
+ Mhat <<- c(M1, M2, M3)
n <<- nrow(X)
d <<- length(M1)
computeW = function(θ)
{
- dim <- d + d^2 + d^3
- W <<- solve( matrix( .C("Compute_Omega",
- X=as.double(X), Y=as.double(Y), M=as.double(M(θ)),
- pn=as.integer(n), pd=as.integer(d),
- W=as.double(W), PACKAGE="morpheus")$W, nrow=dim, ncol=dim) )
+ #require(MASS)
+ dd <- d + d^2 + d^3
+ W <<- MASS::ginv( matrix( .C("Compute_Omega",
+ X=as.double(X), Y=Y, M=Moments(θ), pn=as.integer(n), pd=as.integer(d),
+ W=as.double(W), PACKAGE="morpheus")$W, nrow=dd, ncol=dd ) )
NULL #avoid returning W
},
- M <- function(θ)
+ Moments = function(θ)
{
"Vector of moments, of size d+d^2+d^3"
β2 <- apply(β, 2, function(col) col %o% col)
β3 <- apply(β, 2, function(col) col %o% col %o% col)
- matrix(c(
+ c(
β %*% (p * .G(li,1,λ,b)),
β2 %*% (p * .G(li,2,λ,b)),
- β3 %*% (p * .G(li,3,λ,b))), ncol=1)
+ β3 %*% (p * .G(li,3,λ,b)))
},
f = function(θ)
{
- "Product t(Mi - hat_Mi) W (Mi - hat_Mi) with Mi(theta)"
+ "Product t(hat_Mi - Mi) W (hat_Mi - Mi) with Mi(theta)"
- A <- M(θ) - Mhat
+ L <- expArgs(θ)
+ A <- as.matrix(Mhat - Moments(L))
t(A) %*% W %*% A
},
{
"Gradient of f, dimension (K-1) + d*K + K = (d+2)*K - 1"
- -2 * t(grad_M(θ)) %*% getW(θ) %*% (Mhat - M(θ))
- }
+ L <- expArgs(θ)
+ -2 * t(grad_M(L)) %*% W %*% as.matrix((Mhat - Moments(L)))
+ },
grad_M = function(θ)
{
"Gradient of the vector of moments, size (dim=)d+d^2+d^3 x K-1+K+d*K"
- L <- expArgs(θ)
- p <- L$p
- β <- L$β
+ p <- θ$p
+ β <- θ$β
λ <- sqrt(colSums(β^2))
μ <- sweep(β, 2, λ, '/')
- b <- L$b
+ b <- θ$b
res <- matrix(nrow=nrow(W), ncol=0)
# Gradient on p: K-1 columns, dim rows
km1 = 1:(K-1)
res <- cbind(res, rbind(
- t( sweep(as.matrix(β [,km1]), 2, G1[km1], '*') - G1[K] * β [,K] ),
- t( sweep(as.matrix(β2[,km1]), 2, G2[km1], '*') - G2[K] * β2[,K] ),
- t( sweep(as.matrix(β3[,km1]), 2, G3[km1], '*') - G3[K] * β3[,K] )))
+ sweep(as.matrix(β [,km1]), 2, G1[km1], '*') - G1[K] * β [,K],
+ sweep(as.matrix(β2[,km1]), 2, G2[km1], '*') - G2[K] * β2[,K],
+ sweep(as.matrix(β3[,km1]), 2, G3[km1], '*') - G3[K] * β3[,K] ))
- # TODO: understand derivatives order and match the one in optim init param
for (i in 1:d)
{
# i determines the derivated matrix dβ[2,3]
dβ3_right[block,] <- dβ3_right[block,] + β2
dβ3 <- dβ3_left + sweep(dβ3_right, 2, p * G3, '*')
- res <- cbind(res, rbind(t(dβ), t(dβ2), t(dβ3)))
+ res <- cbind(res, rbind(dβ, dβ2, dβ3))
}
# Gradient on b
res <- cbind(res, rbind(
- t( sweep(β, 2, p * G2, '*') ),
- t( sweep(β2, 2, p * G3, '*') ),
- t( sweep(β3, 2, p * G4, '*') )))
+ sweep(β, 2, p * G2, '*'),
+ sweep(β2, 2, p * G3, '*'),
+ sweep(β3, 2, p * G4, '*') ))
res
},
else if (any(is.na(θ0$b)))
stop("θ0$b cannot have missing values")
- op_res = constrOptim( linArgs(θ0), .self$f, .self$grad_f,
- ui=cbind(
- rbind( rep(-1,K-1), diag(K-1) ),
- matrix(0, nrow=K, ncol=(d+1)*K) ),
- ci=c(-1,rep(0,K-1)) )
-
- # debug:
- print(computeW(expArgs(op_res$par)))
- # We get a first non-trivial estimation of W
- # TODO: loop, this redefine f, so that we can call constrOptim again...
- # Stopping condition? N iterations? Delta <= ε ?
+ # TODO: stopping condition? N iterations? Delta <= epsilon ?
+ for (loop in 1:10)
+ {
+ op_res = constrOptim( linArgs(θ0), .self$f, .self$grad_f,
+ ui=cbind(
+ rbind( rep(-1,K-1), diag(K-1) ),
+ matrix(0, nrow=K, ncol=(d+1)*K) ),
+ ci=c(-1,rep(0,K-1)) )
+
+ computeW(expArgs(op_res$par))
+ # debug:
+ #print(W)
+ print(op_res$value)
+ print(expArgs(op_res$par))
+ }
expArgs(op_res$par)
}
projects / morpheus.git / blobdiff