I am working on the problem of numerically integrating a function $q(\mathbf{x})$ for $\mathbf{x} \in \mathbb{R}^m$, $q:\mathbb{R}^m \rightarrow \mathbb{R}$. The exact form of $q(\mathbf{x})$ will never be known a priori however is, to the first-order, very closely approximated by $$ q(\mathbf{x}) \approx \sum_{i=1}^N c(\mathbf{x};\mathbf{\rho}_{i},\sigma_i,m) \label{eq1} $$ where $\mathbf{\rho}_i \in \mathbb{R}^m$, $\gamma \in \mathbb{R}^+$, $\sigma_i:= \sigma(\mathbf{\rho}_i;\gamma) \in \mathbb{R}^+$ are the location and scale parameters, respectively, of the pdf $c(\mathbf{x};\mathbf{\rho},\sigma,m):\mathbb{R}^m \rightarrow \mathbb{R}$ of an arbitrary multivariate probability distribution. Whilst $\gamma>0$ is a known (small) constant, the values of $\mathbf{\rho}_i$ are unkown a priori.
The problem therefore is to compute numerically $$ I = \int_{a_m}^{b_m} \cdots \int_{a_1}^{b_1} q(\mathbf{x}) \, \text{d}x_1...\text{d}x_m, $$ for finite $(a_j,b_j)$, $j=1,...,m$, as efficiently as possible given the above knowledge of the approximate form of $q(\mathbf{x})$. In particular, that
(i) since $c(\mathbf{x};\mathbf{\rho}_{i},\sigma_i,m)$ is a pdf, much of the computational effort can be expected to be concentrated in the (small) vicinities of the locations $\mathbf{\rho}_{i}$, $i=1,...N$, with lesser effort expended elsewhere,
(ii) we have that $\int_{-\infty}^{+\infty} c(\mathbf{x};\mathbf{\rho}_i,\sigma_i,m) \, \text{d}\mathbf{x} = 1$
(iii) $\arg \max_{\mathbf{x} \in \mathbb{R}^m} c(\mathbf{x};\mathbf{\rho}_{i},\sigma_i,m) = \mathbf{\rho}_i$ at which point $c(\mathbf{\rho}_i;\mathbf{\rho}_{i},\sigma_i,m) = c_{\text{max}}(\sigma_i,m)$, therefore $0 \leq \max_{\mathbf{x} \in \mathbb{R}^m} q(\mathbf{x}) \leq \sum_{i=1}^{N}c_{\text{max}}(\sigma_i,m)$
Some practical notes due to the applications:
(i) $m$ is not large. The ambition is to achieve a scheme that is implementable perhaps for $m \leq 5$,
(ii) limits of integration $(a_j,b_j)$ for $j=1,...m$ collectively describe a hypercube $\mathcal{D}^m$ rather than any arbitrary region,
(iii) $q(\mathbf{x})$ is smooth, continuous and differentiable within $\mathcal{D}^m$,
(iv) it can be assumed that the multiple integral may be evaluated as an iterated integral,
(v) $\gamma$ should be set as small as is feasibly possible as to not hinder the computation, mindful of the effect of concentrating the functional variation close to each $\mathbf{\rho}_i$,
(vi) though theoretically unbounded, $N$ can be assumed to be finite and small (consider $N \approx 2^m$).
Because the numerical evaluation of $q(\mathbf{x})$ may be arbitrarily computationally expensive I therefore wish to ask the community what opinions it has on either or both of the following:
- Which quadrature rule is best suited to efficiently solving a problem of this kind given the knowledge of the approximate functional form of $q(\mathbf{x})$?
- What transformations or scaling of the function $q(\mathbf{x})$ may result in fewer iterations if an adaptive quadrature scheme is employed?
I am grateful for all engagement and of course, if there are any ambiguities in the phrasing of the question, please advise so I can correct for the benefit of all participants.