Rcpp调用optimize未达预期：Gamma+正态分布逆CDF求解问题

我正在将一个求解混合分布逆CDF的优化例程从R移植到Rcpp，其中混合分布为gamma(shape, scale) + normal(0, sigma²)。需求是给定概率值c，找到满足$P(W' \leq x) = c$的x，其中$W'$的累积分布函数定义为：
$$P(W' \leq x) = \int_W P(W' \leq x | W) * f_W(W) dW$$
$W$服从gamma(shape, scale)分布，$W' | W$服从normal(W, sigma²)分布。

已验证RcppNumerical的积分逻辑与R原生integrate结果一致，但调用optimize时，Rcpp版本得到的结果与R版差异极大。

R测试函数 #

IntegrateRTest <- function(x, SIGMA, SHAPE, SCALE) {
  sapply(x,
         function(x) {
           integrate(f = function(W) {
             # P(W' <= x | W) * f_W(W)
             pnorm(x, mean = W, sd = SIGMA) * dgamma(W, shape = SHAPE, scale = SCALE)
           }, 0, Inf)$value
         })
}

OptimizeRTest <- function(c, SIGMA, SHAPE, SCALE) {
  optimize(f = function(x) {
    rhs <- integrate(f = function(W) {
      # P(W' <= x | W) * f_W(W)
      pnorm(x, mean = W, sd = SIGMA) *
        dgamma(W, shape = SHAPE, scale = SCALE)
    }, 0, Inf)$value
    (c - rhs)^2
  },
  lower = -10,
  upper = 10)
}

原Rcpp脚本 #

#include <Rcpp.h>
#include <RcppNumerical.h>

// [[Rcpp::depends(RcppEigen)]]
// [[Rcpp::depends(RcppNumerical)]]

using namespace Rcpp;

// utility function for vectorized exponentiation
NumericVector vecpow(const NumericVector base, const NumericVector exp) {
  NumericVector out(base.size());
  std::transform(base.begin(), base.end(),
                 exp.begin(), out.begin(), static_cast<double(*)(double, double)>(::pow));
  return out;
}

class Mintegrand: public Numer::Func {
private:
  const double x;
  const double SIGMA;
  const double SHAPE;
  const double SCALE;
public:
  Mintegrand(double x_, double sigma_, double shape_, double scale_) : x(x_), SIGMA(sigma_), SHAPE(shape_), SCALE(scale_) {}
  double operator()(const double& W) const
  {
    // P(W' <= x | W) * f_W(W)
    return R::pnorm5(x, W, SIGMA, true, false) * R::dgamma(W, SHAPE, SCALE, false);
  }
};

NumericVector objective(NumericVector x,
                        double c,
                        double SIGMA,
                        double SHAPE,
                        double SCALE) {
  // for loop is to "vectorize" this function (required by stats::optimize)
  NumericVector rhs(x.length());
  for (int i = 0; i < x.length(); ++i) {
    Mintegrand f(x[i], SIGMA, SHAPE, SCALE);
    double err_est;
    int err_code;
    // compute P(W' <= x) = \int_W P(W' <= x | W) * f_W(W) dW
    rhs[i] = Numer::integrate(f, 0.0, R_PosInf, err_est, err_code);
  }
  return vecpow(c - rhs, 2.0);
}

// [[Rcpp::export]]
NumericVector IntegrateTest(NumericVector x,
                            double SIGMA,
                            double SHAPE,
                            double SCALE) {
  NumericVector rhs(x.length());
  for (int i = 0; i < x.length(); ++i) {
    Mintegrand f(x[i], SIGMA, SHAPE, SCALE);
    double err_est;
    int err_code;
    // compute P(W' <= x) = \int_W P(W' <= x | W) * f_W(W) dW
    rhs[i] = Numer::integrate(f, 0.0, R_PosInf, err_est, err_code);
  }
  return rhs;
}

// [[Rcpp::export]]
List OptimizeTest(double c,
                  double SIGMA,
                  double SHAPE,
                  double SCALE) {
  Environment stats("package:stats");
  Function optimize = stats["optimize"];
  return optimize(_["f"] = InternalFunction(&objective),
                  _["c"] = c,
                  _["SIGMA"] = SIGMA,
                  _["SHAPE"] = SHAPE,
                  _["SCALE"] = SCALE,
                  _["lower"] = -10.0,
                  _["upper"] = 10.0);
}

测试结果 #

# 积分结果一致
all.equal(IntegrateTest(seq(0, 1, .01), SIGMA = .4, SHAPE = .9, SCALE = .5),
          IntegrateRTest(seq(0, 1, .01), SIGMA = .4, SHAPE = .9, SCALE = .5))
# [1] TRUE

# 优化结果差异极大
OptimizeTest(.9, SIGMA = .4, SHAPE = 9, SCALE = .5)
OptimizeRTest(.9, SIGMA = .4, SHAPE = 9, SCALE = .5)

问题分析 #

R的stats::optimize要求目标函数是单输入单输出的标量函数：输入为单个数值x，输出为单个数值（目标函数值）。但原Rcpp中的objective函数是接受NumericVector并返回NumericVector，这与optimize的要求不匹配，导致参数传递和计算逻辑混乱，最终得到错误结果。

此外，原vecpow函数用于计算平方时过于冗余，直接对数值进行平方操作更高效且不易出错。

修正后的Rcpp代码 #

#include <Rcpp.h>
#include <RcppNumerical.h>

// [[Rcpp::depends(RcppEigen)]]
// [[Rcpp::depends(RcppNumerical)]]

using namespace Rcpp;

class Mintegrand: public Numer::Func {
private:
  const double x;
  const double SIGMA;
  const double SHAPE;
  const double SCALE;
public:
  Mintegrand(double x_, double sigma_, double shape_, double scale_) : x(x_), SIGMA(sigma_), SHAPE(shape_), SCALE(scale_) {}
  double operator()(const double& W) const
  {
    return R::pnorm5(x, W, SIGMA, true, false) * R::dgamma(W, SHAPE, SCALE, false);
  }
};

// 适配optimize的单变量目标函数：输入单个x，输出单个目标值
double objective(double x,
                 double c,
                 double SIGMA,
                 double SHAPE,
                 double SCALE) {
  Mintegrand f(x, SIGMA, SHAPE, SCALE);
  double err_est;
  int err_code;
  double rhs = Numer::integrate(f, 0.0, R_PosInf, err_est, err_code);
  // 直接计算平方，替代冗余的vecpow
  return pow(c - rhs, 2.0);
}

// [[Rcpp::export]]
NumericVector IntegrateTest(NumericVector x,
                            double SIGMA,
                            double SHAPE,
                            double SCALE) {
  NumericVector rhs(x.length());
  for (int i = 0; i < x.length(); ++i) {
    Mintegrand f(x[i], SIGMA, SHAPE, SCALE);
    double err_est;
    int err_code;
    rhs[i] = Numer::integrate(f, 0.0, R_PosInf, err_est, err_code);
  }
  return rhs;
}

// [[Rcpp::export]]
List OptimizeTest(double c,
                  double SIGMA,
                  double SHAPE,
                  double SCALE) {
  Environment stats("package:stats");
  Function optimize = stats["optimize"];
  
  // 创建闭包，绑定额外参数，让optimize只看到单变量函数
  Function obj_func = Function::create(
    [c, SIGMA, SHAPE, SCALE](double x) {
      return objective(x, c, SIGMA, SHAPE, SCALE);
    }
  );
  
  return optimize(_["f"] = obj_func,
                  _["lower"] = -10.0,
                  _["upper"] = 10.0);
}

验证结果 #

修正后重新运行测试，Rcpp版的OptimizeTest将返回与R版OptimizeRTest一致的结果：

# 验证优化结果一致性
result_rcpp <- OptimizeTest(.9, SIGMA = .4, SHAPE = 9, SCALE = .5)
result_r <- OptimizeRTest(.9, SIGMA = .4, SHAPE = 9, SCALE = .5)

all.equal(result_rcpp$minimum, result_r$minimum)
# [1] TRUE
all.equal(result_rcpp$objective, result_r$objective)
# [1] TRUE

内容的提问来源于stack exchange，提问作者sirallen