Calling optimize in Rcpp not producing expected result

I've been trying to port an optimization routine from R to Rcpp. The Rcpp version isn't producing what I expect and I'm stumped as to what the problem might be. For context, the problem is to compute the inverse cdf of a "gamma(shape, scale) + normal(0, sigma^2)" distribution. In particular, given a value c, find x such that P(W' <= x) = c, where W' has the distribution described. Notice that P(W' <= x) = \int_W P(W' <= x | W) * f_W(W) dW, where W ~ gamma(shape, scale) and W' | W ~ normal(W, sigma^2).

I'm using RcppNumerical (https://stackoverflow.com/a/39449199/2875572) for integration (this seems to be working fine, as the test results indicate). It's the call to optimize that is producing the mysterious results.

R test functions:

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)
}

The Rcpp script:

#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);
}

Test results:

all.equal(IntegrateTest(seq(0, 1, .01), SIGMA = .4, SHAPE = .9, SCALE = .5),
          IntegrateRTest(seq(0, 1, .01), SIGMA = .4, SHAPE = .9, SCALE = .5))
# TRUE

OptimizeTest(.9, SIGMA = .4, SHAPE = 9, SCALE = .5)
OptimizeRTest(.9, SIGMA = .4, SHAPE = 9, SCALE = .5)

# gives very different results