Числовое интегрирование GSL для функции без дополнительных параметров - PullRequest
Новая
       6

Числовое интегрирование GSL для функции без дополнительных параметров

0 голосов
/ 27 января 2020

У меня проблема в том, что мне нужно численно интегрировать одномерную функцию с несколькими дополнительными входами, отличными от переменной, в которую интегрируется. Интеграция - от нуля до бесконечности.

Я сказал без дополнительных параметров , потому что я уже определил класс с дополнительными параметрами, являющимися частными переменными-членами. И тогда определяется функтор оператора, который принимает только переменную интегрирования (следовательно, одномерный). В этом классе я хочу использовать библиотеку числовой интеграции 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 = &alpha; // I don't need this part

        gsl_integration_qagiu (&F, 0.0, 0, 1, 0, 1e-7, 1000,
                        w, &result, &error);
        return result;
    }
};

1 Ответ

0 голосов
/ 28 января 2020

Я нашел решение для этого. Решением было отойти от GSL и использовать библиотеку Boost. В математической библиотеке Boost есть квадратурная функция Гаусса-Кронрода, поэтому она сделает эту работу. 

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