fsbmat fsbmat -4 years ago 178
R Question

Return of the gradient function when using Optim functions and colSums

I can not understand why when using the

gradlik
function as argument to the
Optim
function I get the following error:

Error in optim(beta, loglik, gradlik, method = "BFGS", hessian = T, control = list(fnscale = -1)):
gradient in optim evaluated to length 9000 not 9


However, by calling the
gradlik (beta)
function it returns the gradient vector as expected!

Does anyone have any suggestions for correcting this code?

loglik <- function(beta) {
NXS <- dim(model.matrix(~XS))[2]#Numbers of columns of XS+1
NXO <- dim(model.matrix(~XO))[2]#Numbers of columns of XO+1
## parameter indices
ibetaS <- 1:NXS
ibetaO <- seq(tail(ibetaS, 1)+1, length=NXO)
isigma <- tail(ibetaO, 1) + 1
irho <- tail(isigma, 1) + 1
g <- beta[ibetaS]
b <- beta[ibetaO]
sigma <- beta[isigma]
if(sigma < 0) return(NA)
rho <- beta[irho]
if( ( rho < -1) || ( rho > 1)) return(NA)

XS.g <- model.matrix(~XS) %*% g
XO.b <- model.matrix(~XO) %*% b
u2 <- YO - XO.b
r <- sqrt( 1 - rho^2)
B <- (XS.g + rho/sigma*u2)/r
ll <- ifelse(YS == 0,
(pnorm(-XS.g, log.p=TRUE)),
dnorm(u2/sigma, log = TRUE) - log(sigma) +
(pnorm(B, log.p=TRUE))
)
sum(ll)
}

gradlik <- function(beta) {
NXS <- dim(model.matrix(~XS))[2]
NXO <- dim(model.matrix(~XO))[2]
nObs <- length(YS)
NO <- length(YS[YS > 0])
nParam <- NXS + NXO + 2 #Total of parameters

XS0 <- XS[YS==0,,drop=FALSE]
XS1 <- XS[YS==1,,drop=FALSE]
YO[is.na(YO)] <- 0
YO1 <- YO[YS==1]
XO1 <- XO[YS==1,,drop=FALSE]
N0 <- sum(YS==0)
N1 <- sum(YS==1)

w <- rep(1,N0+N1 )
w0 <- rep(1,N0)
w1 <- rep(1,N1)
NXS <- dim(model.matrix(~XS))[2]
NXO <- dim(model.matrix(~XO))[2]
## parameter indices
ibetaS <- 1:NXS
ibetaO <- seq(tail(ibetaS, 1)+1, length=NXO)
isigma <- tail(ibetaO, 1) + 1
irho <- tail(isigma, 1) + 1

g <- beta[ibetaS]
b <- beta[ibetaO]
sigma <- beta[isigma]
if(sigma < 0) return(matrix(NA, nObs, nParam))
rho <- beta[irho]
if( ( rho < -1) || ( rho > 1)) return(matrix(NA, nObs, nParam))
XS0.g <- as.numeric(model.matrix(~XS0) %*% g)
XS1.g <- as.numeric(model.matrix(~XS1) %*% g)
XO1.b <- as.numeric(model.matrix(~XO1) %*% b)
# u2 <- YO1 - XO1.b
u2 <- YO1 - XO1.b
r <- sqrt( 1 - rho^2)
# B <- (XS1.g + rho/sigma*u2)/r
B <- (XS1.g + rho/sigma*u2)/r
lambdaB <- exp( dnorm( B, log = TRUE ) - pnorm( B, log.p = TRUE ) )
gradient <- matrix(0, nObs, nParam)
gradient[YS == 0, ibetaS] <- - w0 * model.matrix(~XS0) *
exp( dnorm( -XS0.g, log = TRUE ) - pnorm( -XS0.g, log.p = TRUE ) )
gradient[YS == 1, ibetaS] <- w1 * model.matrix(~XS1) * lambdaB/r
gradient[YS == 1, ibetaO] <- w1 * model.matrix(~XO1) * (u2/sigma^2 - lambdaB*rho/sigma/r)
gradient[YS == 1, isigma] <- w1 * ( (u2^2/sigma^3 - lambdaB*rho*u2/sigma^2/r) - 1/sigma )
gradient[YS == 1, irho] <- w1 * (lambdaB*(u2/sigma + rho*XS1.g))/r^3
return(colSums(gradient))
}

n=1000
X1 <- runif(n)
X2 <- runif(n)
XO <- cbind(X1,X2)
X3 <- runif(n)
XS <- cbind(X1,X2,X3)
YS <- sample(c(0,1),n,replace = TRUE)
YO <- sample(100:400,n,replace = TRUE)*YS
beta <- c(1,1,1,1,1,1,1,1,0.5)

#Note that the function below compiles normally:

gradlik(beta)

#But the Optim function does not compile:

theta <-optim(beta,loglik, gradlik, method = "BFGS",hessian = T,control=list(fnscale=-1))
theta$par

Answer Source

Your gradient function needs to give as output a vector with the same size as the number of parameters.

While your final return() is indeed a vector, in your current implementation, there are two other return() in the middle of the code where you still return a matrix.

For instance, when sigma <0 your code returns:

if(sigma < 0) return(matrix(NA, nObs, nParam))

Which is a 9000 x 9 matrix, hence making optim() complain as stated in its error message.

Also when ( rho < -1) || ( rho > 1) your function returns:

if( ( rho < -1) || ( rho > 1)) return(matrix(NA, nObs, nParam))

Which, again, is a 9000 x 9 matrix, resulting in the error.

Therefore, you should start fixing those parts of the code, changing them to return a vector with the same size as the number of parameters.

To see an example of your code returning a matrix, run this:

gradlik(rep(-1, 9))
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
   [1,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
   [2,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
   [3,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
   [4,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
   [5,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
   [6,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
   [7,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
   [8,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
   [9,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [10,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [11,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [12,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [13,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [14,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [15,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [16,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [17,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [18,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [19,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [20,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [21,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [22,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [23,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [24,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [25,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [26,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [27,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [28,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [29,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [30,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [31,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [32,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [33,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [34,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [35,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [36,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [37,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [38,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [39,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [40,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [41,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [42,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [43,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [44,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [45,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [46,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [47,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [48,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [49,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [50,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [51,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [52,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [53,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [54,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [55,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [56,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [57,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [58,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [59,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [60,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [61,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [62,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [63,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [64,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [65,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [66,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [67,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [68,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [69,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [70,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [71,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [72,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [73,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [74,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [75,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [76,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [77,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [78,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [79,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [80,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [81,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [82,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [83,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [84,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [85,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [86,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [87,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [88,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [89,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [90,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [91,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [92,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [93,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [94,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [95,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [96,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [97,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [98,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
  [99,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
 [100,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
 [101,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
 [102,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
 [103,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
 [104,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
 [105,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
 [106,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
 [107,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
 [108,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
 [109,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
 [110,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
 [111,]   NA   NA   NA   NA   NA   NA   NA   NA   NA
 [ reached getOption("max.print") -- omitted 889 rows ]
Recommended from our users: Dynamic Network Monitoring from WhatsUp Gold from IPSwitch. Free Download