【发布时间】：2019-01-17 11:28:17
【问题描述】：
我正在尝试获取基本上相当于逻辑回归模型的 MAP 估计值。我正在使用 optim 函数，它采用对数后验密度及其分析梯度作为参数。我有密度和梯度函数的 R 版本和 Rcpp 版本。我可以使用 R 函数成功估计 MAP，但是 optim 正在进入渐近线并且无法使用 Rcpp 函数收敛到最优值。
我已验证密度函数的 R 版本和密度函数的 Rcpp 版本返回相同的值：
ll_cpp = cpp_posterior_density(THETAi = as.vector(THETA0_LF[i,]),
Yi = as.vector(Y[i,]),
MUi =as.vector(MU_LF[i,]),
invS = invS ,TAU = TAU ,LAMBDA = LAMBDA, J = J ,K = K)
ll_R = lf_posterior_density(THETAi = as.vector(THETA0_LF[i,]),
Yi = as.vector(Y[i,]),
MUi =as.vector(MU_LF[i,]),
invS = invS ,TAU = TAU ,LAMBDA = LAMBDA, J = J ,K = K)
print(paste0(c("R: log posterior: ", ll_R)))
print(paste0(c("cpp: log posterior: ", ll_cpp)))
结果
"R: log posterior: " "15.8951804436067"
"cpp: log posterior: " "15.8951804436067"
我还验证了两个版本的渐变是相等的。
d_cpp = grad(cpp_posterior_density, x = as.vector(THETA0_LF[i,]),
Yi = as.vector(Y[i,]),
MUi =as.vector(MU_LF[i,]),
invS = invS ,TAU = TAU ,LAMBDA = LAMBDA, J = J ,K = K)
d_R = grad(lf_posterior_density, x = as.vector(THETA0_LF[i,]),
Yi = as.vector(Y[i,]),
MUi =as.vector(MU_LF[i,]),
invS = invS ,TAU = TAU ,LAMBDA = LAMBDA, J = J ,K = K)
print(paste0(c("R: gradient of log posterior: ", paste(d_R, collapse = ", "))))
print(paste0(c("cpp: gradient of log posterior: ", paste(d_cpp, collapse = ", "))))
结果
[1] "R: gradient of log posterior: "
[2] "6.49720418347811, 4.67847452089852, 5.93682469664212, 1.47670777676947"
[1] "cpp: gradient of log posterior: "
"6.49720418347811, 4.67847452089852, 5.93682469664212, 1.47670777659075"
但是，当我使用 Rcpp 函数调用 optim 时，我无法收敛：
#Using Rcpp
out_LF = optim(par = as.vector(THETA0_LF[i,]),
fn = cpp_posterior_density,
gr = cpp_grad_posterior_density,
Yi = as.vector(Y[i,]),
MUi = as.vector(MU_LF[i,]),
invS =invS,
TAU = TAU,
LAMBDA = LAMBDA,
J = J,
K = K,
method = "BFGS",
hessian = TRUE,
control = list(trace = 6)) #does not converge
结果
initial value 15.895180
final value -4748.586405
最终值必须严格大于零，表示不收敛。但是，使用 R 函数，我确实得到了收敛：
#With R functions for density and gradient
out_LF2 = optim(par = as.vector(THETA0_LF[i,]),
fn = lf_posterior_density,
gr = lf_grad_posterior_density,
Yi = as.vector(Y[i,]),
MUi = as.vector(MU_LF[i,]),
invS =invS,
TAU = TAU,
LAMBDA = LAMBDA,
J = J,
K = K,
method = "BFGS",
hessian = TRUE,
control = list(trace = 6)) #converged
结果
initial value 15.895180
final value 11.980282
有什么线索吗？
为了重现性，这里是a link to a Dropbox folder，其中包含所需的数据（例如，THETA0_LF、Y、MU_LF 等）以及目标函数和梯度（R 版本和 Rcpp 版本）。还包括一个再现上述输出的 R 文件（参见“debug-rcpp-for-credi.R”）。
下面是目标函数的Rcpp版本
#include <RcppArmadillo.h>
using namespace Rcpp;
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
double cpp_posterior_density(const arma::vec& THETAi, const arma::vec& Yi, const arma::vec& MUi, const arma::mat& invS, const arma::vec& TAU, const arma::mat& LAMBDA, const int J, const int K) {
int j;
double lodd_j;
double b;
// PYi
arma::vec LT = LAMBDA*THETAi;
arma::vec PYi(J);
for (j = 0; j < J; j++){
lodd_j = LT(j) - TAU(j);
if(lodd_j<0){
b = 0;
} else {
b = lodd_j;
}
PYi(j) = exp(lodd_j-b)/(exp(-b) + exp(lodd_j-b));
}
double ll = 0.0;
for (j = 0; j < J; j++){
if (Yi(j)==1L){
ll += log(PYi(j));
}
if (Yi(j)==0L){
ll += log(1.0-PYi(j));
}
}
//Prior distriubtion
arma::vec dMUi = THETAi-MUi;
double twoprior = as_scalar(dMUi.t()*invS*dMUi);
// Return result
double dpost = -1.0*ll - 0.5*twoprior;
return dpost;
}
以下是目标函数的R版本：
lf_posterior_density<-function(THETAi, Yi, MUi, invS, TAU, LAMBDA,J,K, weight = NULL){
if (is.null(weight)){weight = rep(1,J)}
# Defined variables
# PYi - J (vector)
# ll - (scalar)
# dMUi - K (vector)
# prior - (scalar)
# Computations
PYi = as.vector(1/(1 + exp(TAU - LAMBDA%*%THETAi))) # J (vector)
# likelihood component
ll = as.numeric(0) #(scalar)
for (j in 1:J){
if (Yi[j] == 1L){ll = ll + weight[j]*log(PYi[j])}
if (Yi[j] == 0L){ll = ll + weight[j]*log(1.0-PYi[j])}
}
# prior distribution component
dMUi = (THETAi - MUi) # K (vector)
prior = as.numeric(-0.5*(dMUi%*%invS%*%dMUi)) #(scalar)
# Return
return(-ll - prior)
}
【问题讨论】：
-
这是不可重现的。提供目标函数的代码。
-
谢谢，@coatless。我已修改我的提交以使其可重现。