GSL 蒙特卡罗积分的被积函数内的任意多精度答案

【问题标题】：Arbitrary multiprecision inside the integrand of a GSL monte carlo integrationGSL 蒙特卡罗积分的被积函数内的任意多精度
【发布时间】：2022-06-18 01:47:04
【问题描述】：

我的目标是将 GSL 蒙特卡罗积分用于使用任意多精度库 (Boost) 的被积函数。我决定使用任意多精度库，因为积分难以达到收敛。这是描述我要编码的实际数学公式。我的猜测是我没有达到收敛，因此NaN 因为和可以得到非常小的值，接近零。

这是我的代码：

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

以下几行定义了我使用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;
  }

当我尝试使用x0 = 0 运行integrate_integrate 时； 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;我得到NaN。为什么会这样？我没想到会出现这个结果，因为我使用 Log-Sum-Exp 来消除可能的下溢。

【问题讨论】：

如有必要，我可以添加我的test.cpp 文件的标题。
你可以使用boost的蒙特卡洛积分：github.com/boostorg/math/blob/develop/include/boost/math/…
虽然大约需要宇宙的年龄才能收敛。 . .
你为什么这么说？我的意思是，宇宙的年龄..
因为收敛是1/sqrt(n)，其中n是函数调用的次数。所以假设你想恢复双精度（ε=10^-16）。那么你需要n>=10^32。假设你的函数调用每个大约需要 100ns，那将需要 10^25 秒。

标签： c++ numerical-integration gsl arbitrary-precision boost-multiprecision