8.1.1. Introduction

In order to compute the following integral numerically

\[H = \int\limits_{a}^{b} f(x)\,dx\]

one represents it as a sum:

\[H \approx \sum\limits_{i=1}^K \omega_i f(x_i),\]

where \(K\) is the order of integration, \(x_i\) is a set of points of size \(K\) on which the function is computed and \(\omega_i\) are the weights.

The formula is valid regardless of the method is used. In case the rectangular integration (left) is used with \(K=3\), then \(x=\{a, a+\Delta/3, a+2\Delta/3\}\) and \(\omega_i=\Delta/3\), where \(\Delta=b-a\). For \(K=3\) and trapezoidal integration \(x=\{a, (a+b)/2, b\}\) and \(\omega=\{\Delta/4, \Delta/2, \Delta/4\}\).

Other methods may be used, such as Gauss-Legendre or Gauss-Kronrod quadratures. The important point here is that regardless of the method used for the given the order \(K\) and limits \(a\) and \(b\) the sets of \(x_i\) and \(\omega_i\) may be defined once.

The typical problem in calculations is when one needs to compute a set of integrals for each bin to produce a histogram:

\[H_k = \int\limits_{a_k}^{a_{k+1}} f(x)\,dx.\]

this integral is approximated in a similar way:

(1)\[\begin{split}H_k \approx \sum\limits_{i_k=1}^{K_k} \omega_{i_k} f(x_{i_k}), \\ x_{i_k} \in [a_i, a_{i+1}].\end{split}\]

Note, that each bin \(k\) may have different integration precision, different number of integration points and thus has its own index \(i_k\).

The rule may be extended to the higher orders. For example, 2-dimensional integral yielding 2d histogram reads as follows:

(2)\[\begin{split}H_{km} &= \int\limits_{a_k}^{a_{k+1}} dx \int\limits_{b_m}^{b_{m+1}} f(x, y)\,dy \approx \sum\limits_{i_k=1}^{K_k} \sum\limits_{j_m=1}^{M_m} \omega_{i_k}^x \omega_{j_m}^y f(x_{i_k}, y_{j_m}), \\ x_{i_k} &\in [a_i, a_{i+1}], \\ y_{j_m} &\in [b_j, b_{j+1}].\end{split}\]

where we have added a set of bins over variable \(y\) with edges \(b_m\) and integration orders \(M_m\).

Given the formulas (1) and (2) we may now define the integration procedure that involves 2+ transformations:

  1. Sampler transformation. This transformation is provided by the bin edges and required integration orders for each bin. The transformation produces the arrays with integration points \(x\) and weights \(\omega\), needed to compute an integral for each bin. The values of \(x\) and \(\omega\) depend on the quadrature method used.

  2. The integrand transformation or calculation chain that implements \(y_i=f(x_i)\). It receives array \(x\) as input. The integrand is provided by the user.

  3. Integrator transformation. The transformations receives the output of the integrand transformation \(y_i\) as well as the integration weights from the sampler. Integrator implements the convolution and produces the histogram \(H\) as the output. Integrator implementation does not depend on actual sampler used.