r - 从 optim 调用 Rcpp 函数
问题描述
我正在尝试获取基本上相当于逻辑回归模型的 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
有什么线索吗?
为了重现性,这里有一个指向 Dropbox 文件夹的链接,其中包含所需的数据(例如,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)
}
解决方案
您的目标函数有所不同:
ll_cpp = cpp_posterior_density(THETAi = 2*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 = 2*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)))
#> [1] "R: log posterior: " "22.495400131601"
print(paste0(c("cpp: log posterior: ", ll_cpp)))
#> [1] "cpp: log posterior: " "16.7463952181814"
我没有调试你的源代码来发现错误。
REPORT = 1
在这种情况下,添加到control
列表中很有用。对于 R,这给出:
initial value 45.707620
iter 2 value 28.881100
iter 3 value 22.426070
iter 4 value 20.145499
iter 5 value 19.922129
iter 6 value 19.805083
iter 7 value 19.684769
iter 8 value 19.684366
iter 9 value 19.684345
iter 10 value 19.684343
iter 10 value 19.684343
final value 19.684343
converged
对于 Rcpp:
initial value 45.707620
iter 2 value 23.059207
iter 3 value -33.279972
iter 4 value -77.878965
iter 4 value -77.878965
iter 5 value -93.872445
iter 5 value -93.872445
iter 6 value -2830.594586
iter 6 value -2830.594586
iter 6 value -2830.594586
final value -2830.594586
converged
推荐阅读
- reactjs - 如何在 React MaterialUI Select 中重新渲染时自动选择默认选项
- react-native - 如何仅在滑动时禁用视图滚动?
- c# - C#如何根据值重新排序列表元组
- ruby-on-rails - Ruby 3.0.1,Rails 6 puma 服务器连接未建立
- python - 从python中的字符串集中删除不需要的字符
- html - 使用 Iron PDF 将 HTML 转换为 PDF 时没有出现边框
- c# - 自定义输出 JSON 中的 TypeNameHandling 值名称
- c# - EventstoreDb -- AppendToStreamAsync “挂起”。为什么?
- caching - 缓存,你如何处理陈旧的数据
- python - 根据分组数据框中组的前两个值获取数据框