I'm using VEGAS integration, specifically the GSL implementation, for some QCD calculations, and I've been investigating the behavior of the algorithm for various numbers of iterations in an attempt to get more accurate results. The value produced by the integration routine changes rather drastically with the number of iterations, as shown in this plot:
The results shown are for one specific random seed, but they are indicative of other results obtained by running many instances of the program with different random seeds and averaging; specifically, the value with a million iterations is significantly larger than that with a hundred thousand iterations, as shown in this plot (which includes the above value, plus a few other terms whose values are nearly constant):
My best guess as to the cause is that there are peaks in the function being integrated which have a characteristic size around $V/10^6$, with $V$ being the integration volume (in this case, an eight-dimensional hypercube). These peaks would then be detected when sampling with a million points, but sampling with only a hundred thousand points would have a good chance of missing them. That would explain why some of the values at $10^5$ iterations are larger than the values with $10^6$ iterations.
Preliminary question: am I on the right track with this explanation?
Assuming I am, that brings up the question of how to know when I can safely stop increasing the number of iterations. For example, is there some sort of analysis I can apply to the analytic form of the function to pick out the scale of features which will significantly affect the result of the integration?
The integral in question, in case it helps:
$$\begin{multline} \int_{\tau}^1\mathrm{d}z\int_{\tau/z}^1\mathrm{d}\xi\iiint\mathrm{d}^2\vec{q}_1\mathrm{d}^2\vec{q}_2\mathrm{d}^2\vec{q}_3 \\ \frac{1}{1 - \xi}\biggl[\frac{\frac{\tau}{z\xi}g\bigl(\frac{\tau}{z\xi}\bigr)D(z)}{z^2}\frac{(1 - \xi(1-\xi))^2}{\xi^2}F\biggl(\frac{\vec{k}}{\xi} + \vec{q}_3 - \vec{q}_2, z\biggr) - \frac{\frac{\tau}{z}g\bigl(\frac{\tau}{z}\bigr)D(z)}{z^2}F\bigl(\vec{k} + \vec{q}_3 - \vec{q}_2, z\bigr)\biggr] \\ \times F\bigl(\vec{k} + \vec{q}_3 - \vec{q}_1, z\bigr)F(\vec{q}_3, z)\frac{\vec{q}_1\cdot\vec{q}_2}{q_1^2q_2^2} \end{multline}$$ where $F(\vec{q}, z) = \frac{z^{0.288}}{C\pi}\exp\bigl(-q^2 z^{0.288}/C\bigr)$ with $C$ a constant.
