0
$\begingroup$

I'm having a huge problem with Monte Carlo integration. This simple example adequately shows the problem.

ParallelTable[
  NIntegrate[
    3 Sin[θ2] Sqrt[(E^w2/((1 - E^w2 Cos[θ2])^2 + E^(2 w2) Sin[θ2]^2))]
      Sin[θ3] Sqrt[(E^w3/((1 - E^w3 Cos[θ3])^2 + E^(2 w3) Sin[θ3]^2))], 
    {w2, -W, W}, {θ2, 0, Pi}, {w3, -W, W}, {θ3, 0, Pi}, 
    Method -> "QuasiMonteCarlo"] /
  (NIntegrate[
     Sin[θ2]Sqrt[(E^w2/((1 - E^w2 Cos[θ2])^2 + E^(2 w2) Sin[θ2]^2))], 
     {w2, -W, W}, {θ2, 0, Pi},
     Method -> "QuasiMonteCarlo"])^2, 
  {W, 1, 10000, 100}]

This gives

{3.0014, 2.9491, 2.85273, 2.94142, 3.13837, 3.3748, 3.59702, 3.76732, 3.866, 
 3.88818, 3.83872, 3.72786, 3.56807, 3.37201, 3.15141, 2.91649, 2.67575, 2.43598, 
 2.20243, 1.97891, 1.76811, 1.5717, 1.39061, 1.22513, 1.07511, 0.940038, 0.819179, 
 0.711636, 0.616418, 0.532496, 0.458835, 0.394427, 0.338305, 0.28956, 0.247351, 
 0.2109, 0.179505, 0.152529, 0.129403, 0.109619, 0.0927273, 0.0783328, 0.0660878, 
 0.0556886, 0.046871, 0.0394056, 0.033094, 0.027765, 0.0232714, 0.0194869, 
 0.0163033, 0.0136281, 0.0113825, 0.00949942, 0.00792187, 0.00660148, 0.00549731, 
 0.00457473, 0.00380451, 0.00316199, 0.00262639, 0.00218025, 0.00180887, 0.00149994, 
 0.00124312, 0.00102975, 0.000852585, 0.000705567, 0.000583633, 0.000482557, 
 0.000398814, 0.000329467, 0.000272068, 0.000224581, 0.000185313, 0.000152854, 
 0.000126035, 0.000103886, 0.0000856005, 0.0000705103, 0.0000580618, 0.0000477963, 
 0.0000393339, 0.0000323603, 0.0000266154, .0000218844, 0.0000179894, 0.0000147838, 
 0.0000121462, 9.97676*10^-6, 8.19278*10^-6, 6.7262*10^-6, 5.52086*10^-6, 
 4.53049*10^-6, 3.71696*10^-6, 3.04885*10^-6, 2.5003*10^-6, 2.05002*10^-6, 
 1.6805*10^-6, 1.37731*10^-6}

The integral in the numerator is just the integral inside the denominator squared. So this expression is essentially computing $3\,x\,/\,x$. Monte Carlo integration should be calculating the numerator integral to have the same value as the denominator integral squared regardless of how large W is. Thus, you would expect the output to be a list of numbers near 3 if not exactly 3.

But as the results below show, as W increases the Monte Carlo integration becomes worse and worse at estimating the numerator. It progressively underestimates the value of the numerator as W becomes larger and larger. If anyone can help me successfully integrate these types of functions for large W, I would greatly appreciate it.

What I describe here happens even when I use regular Monte Carlo integration and not just Quasi-Monte Carlo.

$\endgroup$
  • 1
    $\begingroup$ I think that the problem is that your integrals only have a significant value over a very small part of their range (that is, when w2 is small). The Monte Carlo picks are probably missing this region. Hitting the region is presumably more difficult in the four-fold integral than in the two-fold integral. $\endgroup$ – mikado Sep 26 '18 at 20:02
  • $\begingroup$ I suspect that as well, they are falling off to quickly. Is there anyway I can amend this? $\endgroup$ – Daniel Berkowitz Sep 26 '18 at 20:04
  • 1
    $\begingroup$ You might try a change of variable, e.g. u==Exp[w]. For accurate Monte Carlo integration, I think you want to try to minimise the variance of the integrand. $\endgroup$ – mikado Sep 26 '18 at 21:12
  • $\begingroup$ I'll try thanks and let you know how it goes on this thread. $\endgroup$ – Daniel Berkowitz Sep 26 '18 at 21:19
  • 2
    $\begingroup$ You can always try the AdaptiveQuasiMonteCarlo Method. This hopefully spends more of the evaluations near the peaks of the function. For more ideas the reference has a great guide to numerical integration with many good examples and also a big section how to tune MonteCarlo integration. $\endgroup$ – Thies Heidecke Sep 27 '18 at 3:32
2
$\begingroup$

If you slightly modify the NIntegrate range specifications of your loop variable W then you'll get results you expect:

AbsoluteTiming[
 Block[{pg = 4}, 
  res = Table[
    NIntegrate[
      3 Sin[\[Theta]2] Sqrt[(E^
           w2/((1 - E^w2 Cos[\[Theta]2])^2 + 
            E^(2 w2) Sin[\[Theta]2]^2))] Sin[\[Theta]3] Sqrt[(E^
           w3/((1 - E^w3 Cos[\[Theta]3])^2 + 
            E^(2 w3) Sin[\[Theta]3]^2))], {w2, -W, -10, 10, 
       W}, {\[Theta]2, 0, Pi}, {w3, -W, -10, 10, W}, {\[Theta]3, 0, 
       Pi}, Method -> "MonteCarlo", 
      PrecisionGoal -> 
       pg]/(NIntegrate[
        Sin[\[Theta]2] Sqrt[(E^
             w2/((1 - E^w2 Cos[\[Theta]2])^2 + 
              E^(2 w2) Sin[\[Theta]2]^2))], {w2, -W, -10, 10, 
         W}, {\[Theta]2, 0, Pi}, Method -> "MonteCarlo", 
        PrecisionGoal -> pg])^2, {W, 1, 10000, 100}]
  ]
 ]

(* {15.1931, {2.19657, 3.05144, 2.99393, 3.03217, 3.01784, 
  3.08816, 3.14745, 2.91172, 3.1147, 3.02755, 2.92361, 2.98091, 
  3.18165, 3.06744, 3.0053, 3.02033, 3.03638, 2.97478, 2.91792, 
  2.89669, 2.99929, 3.09445, 3.05545, 3.03786, 2.97139, 3.07607, 
  3.00059, 3.09282, 2.90901, 3.04438, 3.00764, 3.03369, 2.98443, 
  2.93917, 2.96088, 3.03941, 3.06992, 2.9463, 3.03358, 2.88547, 
  3.04671, 2.89731, 2.91739, 2.99103, 3.05259, 3.15771, 2.87781, 
  3.02416, 3.09235, 2.91561, 2.92959, 3.03615, 3.09252, 2.90106, 
  3.06102, 3.01301, 3.05854, 2.9825, 3.01436, 2.95519, 2.88678, 
  3.15286, 3.01529, 2.92353, 3.04329, 2.98486, 3.05863, 2.94734, 
  3.05641, 2.88573, 2.99265, 3.09635, 3.02641, 2.98033, 2.93819, 
  2.96236, 3.04147, 2.90924, 2.90145, 3.00857, 2.98786, 3.16616, 
  2.97454, 3.05536, 2.99765, 2.92687, 3.01301, 2.94392, 2.95414, 
  3.05081, 2.97163, 2.98609, 2.86618, 2.96695, 3.04867, 3.04276, 
  2.98415, 3.05364, 3.07536, 3.04569}}  *)
$\endgroup$
  • $\begingroup$ Thanks a lot for the help, sorry for not replying sooner, I'm working with these types of calculations in my dissertation for my Ph.D in theoretical physics, so I was really absorbed in my work for the past couple of days. The actual integral I'm calculating falls off much faster, the method you outlined has helped but I'm still getting major discrepancies. So I have a better understanding of how NIntegrate works by adding the -10 and 10 to the range what does that do exactly? $\endgroup$ – Daniel Berkowitz Sep 29 '18 at 5:05
  • $\begingroup$ NIntegrate's integration range specifications can take intermediate points. The integration region is partitioned at those intermediate points. MinRecursion does that kind of partitioning too, but it is using a regular grid specified by a single integer. $\endgroup$ – Anton Antonov Sep 29 '18 at 13:58

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.