0

My aim is to use GSL monte carlo integration for an integrand in which an arbitrary multiprecision library (Boost) is used. I decided to use an arbitrary multiprecision library because the integral had difficulties to reach convergence. This is the actual mathematical formula describing what I am trying to code. My guess is that I do not reach convergence and therefore NaN because and can get very small values, near zeros.

This is my code:

mp::float128 PDFfunction(double invL, int t, double invtau, double x0, double x, int n_lim) {
    
    const double c =  M_PI * (M_PI/4) * ((2 * t) * invtau);
    mp::float128 res = 0;
    
    for(int n = 1; n <= n_lim; ++n){
      res += exp(-1 * (n * n) * c) * cos((n * M_PI * x) * invL) * cos((n * M_PI * x0) * invL);
    }
    mp::float128 res_tot = invL + ((2 * invL) * res);
return res_tot;
        }

The following lines define the integral that I carry out using GSL:

struct my_f_params {double x0; double xt_pos; double y0; double yt_pos; double invLx; double invLy;
    double invtau_x; double invtau_y; int n_lim; double tax_rate;};


double g(double *k, size_t dim, void *p){

  struct my_f_params * fp = (struct my_f_params *)p;

  mp::float128 temp_pbx = prob1Dbox(fp->invLx, k[0], fp->invtau_x, fp->x0, fp->xt_pos, fp->n_lim);
  mp::float128 temp_pby = prob1Dbox(fp->invLy, k[0], fp->invtau_y, fp->y0, fp->yt_pos, fp->n_lim);
  mp::float128 AFac = (-2 * k[0] * fp->tax_rate);
  mp::float128 res = exp(log(temp_pbx) + log(temp_pby) + AFac);
  return res.convert_to<double>();
  }



double integrate_integral(const double& x0, const double& xt_pos, const double& y0,
    const double& yt_pos, const double& invLx, const double& invLy, const double& invtau_x,
    const double& invtau_y, const int& n_lim, const double& tax_rate){


      double res, err;
      double xl[1] = {0};
      double xu[1] = {10000000};

      const gsl_rng_type *T;
      gsl_rng *r;

      gsl_monte_function G;

      struct my_f_params params = {x0, xt_pos, y0, yt_pos, invLx, invLy, invtau_x, invtau_y,  n_lim, tax_rate};
      G.f = &g;
      G.dim = 1;
      G.params = &params;

      size_t calls = 10000;

      gsl_rng_env_setup ();

      T = gsl_rng_default;
      r = gsl_rng_alloc (T);


      gsl_monte_vegas_state *s = gsl_monte_vegas_alloc (1);

      gsl_monte_vegas_integrate (&G, xl, xu, 1, 10000, r, s,
                             &res, &err);

  do
        {
          gsl_monte_vegas_integrate (&G, xl, xu, 1, calls/5, r, s,
                                     &res, &err);
        }
      while (fabs (gsl_monte_vegas_chisq (s) - 1.0) > 0.5);

      gsl_monte_vegas_free (s);
      gsl_rng_free (r);

  return res;
  }

When I try to run integral_integrate with x0 = 0; xt_pos = 0; y0 = 0; yt_pos = 10; invLx = invLy = 0.09090909; invtau_x = invtau_y = 0.000661157; n_lim = 1000; tax_rate = 7e-8; I get NaN. Why is this the case? I was not expecting this result since I used Log-Sum-Exp to get rid of possible underflow.

CafféSospeso
  • 1,101
  • 3
  • 11
  • 28
  • If necessary, I can add the headers of my `test.cpp` file. – CafféSospeso Jun 06 '22 at 23:07
  • You can use boost's Monte-Carlo integration: https://github.com/boostorg/math/blob/develop/include/boost/math/quadrature/naive_monte_carlo.hpp – user14717 Jun 07 '22 at 13:49
  • Though it should take roughly the age of the universe to converge . . . – user14717 Jun 07 '22 at 13:50
  • Why you say that? I mean, the age of the universe.. – CafféSospeso Jun 07 '22 at 21:29
  • Because the convergence is 1/sqrt(n), where n is the number of function calls. So suppose you want to recover double precision (ε=10^-16). Then you need n>=10^32. Assuming your function calls take roughly 100ns each, that would take 10^25 seconds. – user14717 Jun 08 '22 at 00:36
  • What would you suggest? may be set a conditional statements in the case of x << 1? – CafféSospeso Jun 16 '22 at 22:44
  • That won't work either. I can't see what you're trying specifically trying to achieve without it being rendered in LaTeX. – user14717 Jun 17 '22 at 14:58
  • Ok, so I'm going to edit my question including the actual equations in LaTeX. Thanks for the suggestion. – CafféSospeso Jun 17 '22 at 17:13
  • Ok, you're in low dimensions so you don't need Monte-Carlo integration. Try sparse grids or just tensor producting a 1D quadrature rule. – user14717 Jun 17 '22 at 18:13
  • Well, in reality my actual code (which I did not show here for simplicity) is given by a four-dimensional integral (reason why I am using a Monte-Carlo integration. I just texted you in the chat: https://chat.stackoverflow.com/rooms/245699/room-for-caffesospeso-and-nick-thompson – CafféSospeso Jun 17 '22 at 18:16
  • 4 dimensions is still low for Monte-Carlo. Sparse grids is your friend. – user14717 Jun 17 '22 at 18:29
  • Let us [continue this discussion in chat](https://chat.stackoverflow.com/rooms/245701/discussion-between-caffesospeso-and-user14717). But are you aware of some popular sparse grid integration library in C++? – CafféSospeso Jun 17 '22 at 18:31

0 Answers0