I am trying to understand the rebinning algorithm of the VEGAS (original publication (preprint from LKlevin) and implementation notes) Monte Carlo integration. I will try to explain first what I think I understood and then pose my questions.

For simplicity, let's assume we have a 1-dimensional function $f(x)$ that is positive over the whole interval $[0,1]$. This interval is seperated into, let's say, $n$ bins. These bins are initially equally-sized. The bin sizes $\Delta x_i$ define a probability density

$$\rho(x) = \begin{cases} 0 \le x < \Delta x_1 : \frac{1}{n \Delta x_1} \\ \vdots \\ \Delta x_{n-1} \le x < \Delta x_n : \frac{1}{n \Delta x_n} \end{cases} .$$

The bin sizes must add up to the length of the interval in order to make $\rho(x)$ properly normalized:

$$\sum_{i=1}^n \Delta x_i = 1 \quad \Rightarrow \quad \int\limits_0^1 \mathrm{d}x \, \rho(x) = 1.$$

Using $N$ randomly chosen numbers $\{ x_i \}$ from the distribution $\rho(x)$ VEGAS computes an approximation $S^{(1)}$ of the integral:

$$S^{(1)} = \frac{1}{N} \sum_{\{x_i\}} \frac{f(x_i)}{\rho(x_i)} \approx \int\limits_0^1 \mathrm{d}x \, f(x)$$

Up to now this is just importance sampling (I am not interested in stratified sampling) using a variably-sized grid. The interesting piece in VEGAS is the rebinning algorithm, i.e. the algorithm which recomputes the bin sizes depending on the function values accumulated in the iteration before:

  • For each bin the squared function values (?) are summed up (in the original publication the absolute values are summed up).
  • There is also a dampening function applied to each value, to "avoid rapid, destablizing changes".
  • After that, each value is smoothed with the neighboring bins. I guess this also adds some stability to the rebinning algorithm for certain functions (but I can not explain why). Let's call the final values $b_i$.
  • The bin sizes are now set such that every new bin contains approximately the average:

$$\overline{b} = \frac{1}{n} \sum_{i=1}^n b_i$$

This algorithm makes the bins grow where the function is "small" and shrink where the function is "big". Small and big are understood in relation to each other, e.g. the function's maxima are considered "big" and everything else would be considered "smaller". Since the probability for a point $x$ to end up in any bin is equal (VEGAS is written to exhibit this behavior), the function is sampled most where it is biggest and thereby the error is reduced.

Why is it written that way? Why does one not tackle the problem more directly, for example by taking the averaged and binned function as a probability density for the next iteration?

  • $\begingroup$ You should add a link to the original publication. Further, what do you mean with "functions that differ largely in neighboring bins"? What is a small function? What is a big function? $\endgroup$
    – vanCompute
    Nov 26, 2012 at 22:09
  • $\begingroup$ Not quite sure what your question is. Is the problem that the absolute values are summed? $\endgroup$
    – LKlevin
    Apr 12, 2015 at 11:31
  • $\begingroup$ @LKlevin: Well, the squared values are summed which differs from what is written in the original publication (sum absolute values) for one. But I guess this is because it yields a better convergence for most integrands in practice. But why is that? Did they just do that on an empirical basis or is there a proper explanation? I am also not sure why this algorithm could be stable at all if you sum the squares instead of the absolute values. I tried to find a better rebinning algorithm but it seems to be impossible to find a stable one for integrands that are not factorizable. Why is it so good? $\endgroup$
    – cschwan
    Apr 14, 2015 at 9:20

1 Answer 1


Disclaimer: I'm very familiar with Monte Carlo algorithms in general, but not with the VEGAS algorithm specifically.

The VEGAS algorithm is based on the fact that if we knew the exact probability density that minimizes the variance of our Monte Carlo integration, we could get our answer using exactly 1 evaluation of the function f.

This is given by

$$p(\mathbf{x}) = \frac {|f(\mathbf{x})|} {\int_\Omega |f(\mathbf{x})| d\mathbf{x}}$$

We don't know the probability density, and estimating it exactly is not feasible for a high-dimensional function as it would quickly take up a massive amount of memory.

Instead the VEGAS algorithm approximates it as an M-step piecewise constant function.

I'm guessing you don't have access to the full article? The original article does NOT use the squared function, but the absolute one (a pre-print version can be found here).

I hope this helps answer your question. The original article answers most of this, so it might be worthwhile getting


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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