我有一个问题,我需要对一个单变量函数进行数值积分,而不要对要积分的变量进行多个附加输入。积分是从零到无穷大。
我说没有额外的参数,因为我已经定义了一个类,其中额外的参数是私有成员变量。然后将运算符定义为仅接受积分变量(因此,单变量)。对于此类,我想使用GSL数字积分库(gsl / gsl_integration.h)进行积分。有没有一种方法可以使用GSL在类内部定义此集成的成员函数?
#include <cmath>
#include <Rmath.h>
#include <algorithm>
#include <iterator>
#include <RcppArmadillo.h>
#include <progress.hpp>
#include <progress_bar.hpp>
#include <RcppGSL.h>
#include <gsl/gsl_integration.h>
#include <Rdefines.h>
// [[Rcpp::depends(RcppArmadillo, RcppProgress, RcppGSL)]]
using namespace arma;
class ObservedLik
{
private:
const int& Tk;
const arma::vec& resid;
const arma::mat& ZEREZ_S;
const double& nu;
const double& maxll;
public:
ObservedLik(const int& Tk_,
const arma::vec& resid_,
const arma::mat& ZEREZ_S_,
const double& nu_,
const double& maxll_) : Tk(Tk_), resid(resid_), ZEREZ_S(ZEREZ_S_), nu(nu_), maxll(maxll_) {}
double operator()(const double& lam) const {
double loglik = -M_LN_SQRT_2PI * static_cast<double>(Tk) + (0.5 * nu - 1.0) * lam - 0.5 * nu * lam + 0.5 * nu * (std::log(nu) - M_LN2) - R::lgammafn(0.5 * nu);
double logdet_val;
double logdet_sign;
log_det(logdet_val, logdet_sign, ZEREZ_S);
loglik -= 0.5 * (logdet_val + arma::accu(resid % arma::solve(ZEREZ_S, resid)));
/***********************************
subtract by maximum likelihood value
for numerical stability
***********************************/
return std::exp(loglik - maxll);
}
double integrate() {
/* do the integration here */
gsl_integration_workspace * w
= gsl_integration_workspace_alloc (1000);
double result, error;
gsl_function F;
F.function = &f; // make this the operator()
F.params = α // I don't need this part
gsl_integration_qagiu (&F, 0.0, 0, 1, 0, 1e-7, 1000,
w, &result, &error);
return result;
}
};
我找到了解决方案。解决方案是脱离GSL,使用Boost库。 Boost数学库中有一个Gauss-Kronrod正交函数,因此可以完成此工作。
#include <cmath>
#include <Rmath.h>
#include <algorithm>
#include <iterator>
#include <RcppArmadillo.h>
#include <progress.hpp>
#include <progress_bar.hpp>
#include <boost/math/quadrature/gauss_kronrod.hpp>
#include "dic_nmr.h"
// [[Rcpp::depends(RcppArmadillo, RcppProgress, BH)]]
using namespace arma;
using namespace boost::math::quadrature;
class ObservedLik
{
private:
const int& Tk;
const arma::vec& resid;
const arma::mat& ZEREZ_S;
const double& nu;
const double& maxll;
public:
ObservedLik(const int& Tk_,
const arma::vec& resid_,
const arma::mat& ZEREZ_S_,
const double& nu_,
const double& maxll_) : Tk(Tk_), resid(resid_), ZEREZ_S(ZEREZ_S_), nu(nu_), maxll(maxll_) {}
double integrate_()(double lam) const {
double loglik = -M_LN_SQRT_2PI * static_cast<double>(Tk) + (0.5 * nu - 1.0) * lam - 0.5 * nu * lam + 0.5 * nu * (std::log(nu) - M_LN2) - R::lgammafn(0.5 * nu);
double logdet_val;
double logdet_sign;
log_det(logdet_val, logdet_sign, ZEREZ_S);
loglik -= 0.5 * (logdet_val + arma::accu(resid % arma::solve(ZEREZ_S, resid)));
/***********************************
subtract by maximum likelihood value
for numerical stability
***********************************/
return std::exp(loglik - maxll);
}
double integrate() {
/* do the integration here */
double error;
double Q = gauss_kronrod<double, 31>::integrate(integrate_, 0.0, std::numeric_limits<double>::infinity(), 5, 1e-14, &error);
return Q;
}
};