1
在使用，爲什麼當gradlik功能參數傳遞給Optim功能我碰到下面的錯誤我不明白：返回使用Optim功能和colSums
Error in optim(beta, loglik, gradlik, method = "BFGS", hessian = T, control = list(fnscale = -1)):
gradient in optim evaluated to length 9000 not 9
然而，通過調用gradlik (beta)功能它會像預期的那樣返回漸變向量！
有沒有人有任何建議糾正此代碼？
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
我得到一個錯誤與分配給'theta'代碼：....'錯誤的Optim（測試版，loglik ，gradlik，method =「BFGS」，hessian = T，control = list（fnscale = -1））： 漸變優化評估爲長度9000而不是9'。 「編譯漸變」的含義並不清楚。 R是一種解釋型語言。你在做某種未指定的調試嗎？ –
問題是這樣的，梯度大小爲'9000'而不是'9'？問題很明顯！我不做解釋沒有指定！我只是不明白爲什麼這個錯誤！ – fsbmat