Skip to main content
Melhorei a pergunta ao apresentar o erro que recebo e também acrescentando a função gradlik(beta)
Fonte Link
fsbmat
  • 1,4mil
  • 8
  • 22

nãoNão consigo entender porque o gradiente abaixo retorna uma matriz e não um vetorque ao rodarusar a função gradlik como argumento da função Optim. Quando se compila linha a linha recebo o gradienteseguinte erro:

Error in optim(beta, loglik, gradlik, method = "BFGS", hessian = T, control = list(fnscale = -1)) : 
gradiente em optim retorna um objeto de comprimento 9000 ao invés de 9

Porém, obtém-se umao chamar a função gradlik(beta) a mesma me retorna o vetor gradiente como esperado.!

Alguém possui alguma sugestão para a correção deste código?

loglik <- function(beta) {
  NXS <- dim(model.matrix(~XS))[2]#Número de colunas de XS+1
  NXO <- dim(model.matrix(~XO))[2]#Número de colunas de 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]#Número de colunas de XS+1
  NXO <- dim(model.matrix(~XO))[2]#Número de colunas de XO+1
  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]#Número de colunas de XS+1
  NXO <- dim(model.matrix(~XO))[2]#Número de colunas de 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(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 que a função abaixo compila normalmente:
gradlik(beta)
#Porém a função Optim não compila:
theta <-optim(beta,loglik, gradlik, method = "BFGS",hessian = T,control=list(fnscale=-1)) 
theta$par

Note que com os códigos abaixo a função optim compila corretamente! Mas como desejo fazer algumas simulações, quero algo mais geral e não entendo porquê o gradiente retorna uma matriz ao invés de um vetor com o código acima!

library("ssmrob")
data(MEPS2001)
attach(MEPS2001) 
library("sampleSelection")
XO <- cbind(age,female,educ,blhisp,totchr,ins)
XS <- cbind(age,female,educ,blhisp,totchr,ins,income)
YO <- lnambx
YS <- dambexp
fit<-heckit(YS ~ XS, YO ~ XO)
beta <-  coef(fit)[-16]

não consigo entender porque o gradiente abaixo retorna uma matriz e não um vetor ao rodar a função Optim. Quando se compila linha a linha o gradiente, obtém-se um vetor como esperado. Alguém possui alguma sugestão para a correção deste código?

loglik <- function(beta) {
  NXS <- dim(model.matrix(~XS))[2]#Número de colunas de XS+1
  NXO <- dim(model.matrix(~XO))[2]#Número de colunas de 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]#Número de colunas de XS+1
  NXO <- dim(model.matrix(~XO))[2]#Número de colunas de XO+1
  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]#Número de colunas de XS+1
  NXO <- dim(model.matrix(~XO))[2]#Número de colunas de 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(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)

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

Note que com os códigos abaixo a função optim compila corretamente! Mas como desejo fazer algumas simulações, quero algo mais geral e não entendo porquê o gradiente retorna uma matriz ao invés de um vetor com o código acima!

library("ssmrob")
data(MEPS2001)
attach(MEPS2001) 
library("sampleSelection")
XO <- cbind(age,female,educ,blhisp,totchr,ins)
XS <- cbind(age,female,educ,blhisp,totchr,ins,income)
YO <- lnambx
YS <- dambexp
fit<-heckit(YS ~ XS, YO ~ XO)
beta <-  coef(fit)[-16]

Não consigo entender porque que ao usar a função gradlik como argumento da função Optim recebo o seguinte erro:

Error in optim(beta, loglik, gradlik, method = "BFGS", hessian = T, control = list(fnscale = -1)) : 
gradiente em optim retorna um objeto de comprimento 9000 ao invés de 9

Porém, ao chamar a função gradlik(beta) a mesma me retorna o vetor gradiente como esperado!

Alguém possui alguma sugestão para a correção deste código?

loglik <- function(beta) {
  NXS <- dim(model.matrix(~XS))[2]#Número de colunas de XS+1
  NXO <- dim(model.matrix(~XO))[2]#Número de colunas de 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]#Número de colunas de XS+1
  NXO <- dim(model.matrix(~XO))[2]#Número de colunas de XO+1
  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]#Número de colunas de XS+1
  NXO <- dim(model.matrix(~XO))[2]#Número de colunas de 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(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 que a função abaixo compila normalmente:
gradlik(beta)
#Porém a função Optim não compila:
theta <-optim(beta,loglik, gradlik, method = "BFGS",hessian = T,control=list(fnscale=-1)) 
theta$par
Acrescentei um código para melhor entendimento da minha dúvida!
Fonte Link
fsbmat
  • 1,4mil
  • 8
  • 22

Note que com os códigos abaixo a função optim compila corretamente! Mas como desejo fazer algumas simulações, quero algo mais geral e não entendo porquê o gradiente retorna uma matriz ao invés de um vetor com o código acima!

library("ssmrob")
data(MEPS2001)
attach(MEPS2001) 
library("sampleSelection")
XO <- cbind(age,female,educ,blhisp,totchr,ins)
XS <- cbind(age,female,educ,blhisp,totchr,ins,income)
YO <- lnambx
YS <- dambexp
fit<-heckit(YS ~ XS, YO ~ XO)
beta <-  coef(fit)[-16]

Note que com os códigos abaixo a função optim compila corretamente! Mas como desejo fazer algumas simulações, quero algo mais geral e não entendo porquê o gradiente retorna uma matriz ao invés de um vetor com o código acima!

library("ssmrob")
data(MEPS2001)
attach(MEPS2001) 
library("sampleSelection")
XO <- cbind(age,female,educ,blhisp,totchr,ins)
XS <- cbind(age,female,educ,blhisp,totchr,ins,income)
YO <- lnambx
YS <- dambexp
fit<-heckit(YS ~ XS, YO ~ XO)
beta <-  coef(fit)[-16]
Fonte Link
fsbmat
  • 1,4mil
  • 8
  • 22

Retorno da função gradiente ao usar funções Optim e colSums

não consigo entender porque o gradiente abaixo retorna uma matriz e não um vetor ao rodar a função Optim. Quando se compila linha a linha o gradiente, obtém-se um vetor como esperado. Alguém possui alguma sugestão para a correção deste código?

loglik <- function(beta) {
  NXS <- dim(model.matrix(~XS))[2]#Número de colunas de XS+1
  NXO <- dim(model.matrix(~XO))[2]#Número de colunas de 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]#Número de colunas de XS+1
  NXO <- dim(model.matrix(~XO))[2]#Número de colunas de XO+1
  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]#Número de colunas de XS+1
  NXO <- dim(model.matrix(~XO))[2]#Número de colunas de 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(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)

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