Stability of Algorithms in Statistical Modeling

Abstract

In this paper, we investigate algorithms stability for calculation of multidimensional integrals using the statistical modeling methods. We considered issues of the algorithms optimization and we give sufficient conditions for the stability. We apply our approach to both calculation of integral from the regression function and the moments integral calculation. In all our numerical experiences, we used the mt19937 pseudorandom number generator.

1. Introduction

One of the main problem of the statistical modeling method (the Monte Carlo method) is the problem of quality for pseudorandom numbers. In the paper, we consider a task of multidimensional integrals calculation by the statistical modeling method and give sufficient conditions for the stability of this task to quality of pseudorandom numbers. Included results of various numerical experiences with the mt19937 pseudorandom number generator. In our work, we discuss important issues of algorithms optimization in the statistical modeling. In particular, we apply the new approach to the following: a task of finding of integral functionals from solution of boundary-value problems for both the linear [1] or nonlinear [2] elliptic equations (the estimations are given near to a boundary).

The paper is organized as follows: In Section 2, we give the sufficient conditions of stability. Calculation of an integral of very large dimensions is discussed in Section 3. Rare events effect is the subject of Section 4. In Section 5, we describe calculation of integral moments. In Section 6, we apply our approach to calculate an integral of the regression function. In Section 7, we give the conclusion of our studies.

2. Sufficient conditions of stability

Let

$$I={\int }_{D}f\left(x\right)\mathit{dx}$$E1

be the Riemann integral. Here $D$ is a domain of the s-dimensional Euclidean space $R^s$. If the dimension $s$ is large enough then we must use a statistical modeling method. In this case, our integral has form of the mathematical expectation for a random value $\eta =f\left(\xi \right)/p\left(\xi \right)$:

$${\int }_{D}f\left(x\right)\mathit{dx}={\int }_{D}p\left(x\right)\frac{f\left(x\right)}{p\left(x\right)}\mathit{dx}=\mathcal{E}\frac{f\left(x\right)}{p\left(x\right)}=\mathcal{E}\eta .$$E2

Here $p\left(x\right)$ is a density of random variable $\xi \in D$. We put $p\left(x\right)\ne 0$ for $f\left(x\right)\ne 0$, and we say that there exists integral

$${\int }_D\mid f\left(x\right)\mid \mathit{dx}.$$

A variance of the random value $\eta$:

$${\sigma }^2=\mathbf{var}\eta =\mathcal{E}{\eta }^2-{\left(\mathcal{E}\eta \right)}^2={\int }_{D}p\left(x\right){\left[\frac{f\left(x\right)}{p\left(x\right)}\right]}^2\mathit{dx}-I^2={\int }_D\frac{f^2\left(x\right)}{p\left(x\right)}\mathit{dx}-I^2.$$E3

We estimate the mathematical expectation $\mathcal{E}\eta$ by the sum ${\sum }_{i=1}^N{\eta }_i/N$, where ${\eta }_i$ are independent realizations of the random value $\eta$. Suppose ${\sigma }^2$ is finite, and $N$ is large enough; then from the classical central limit theorem, it follows that the random value ${\sum }_{i=1}^N{\eta }_i/N$ has distribution close to the normal distribution with a mathematical expectation $I$, and mean-square deviation $\sigma /\sqrt{N}$. This property above is useful to estimate error, e.g., using the 3$\sigma$ rule. So we have

$$\left|I-\frac{1}{N}\sum _{i=1}^N{\eta }_i\right|\le \frac{3\sigma }{\sqrt{N}}$$E4

with probability 0,997, approximately.

Let us ${\eta }_i$ be realizable sampling values; then the value $\sigma$ is estimated as the following:

$${\sigma }^2\approx \frac{1}{N}\sum _{i=1}^N{\eta }_i^2-{\left(\frac{1}{N}\sum _{i=1}^N{\eta }_i\right)}^2.$$E5

Suppose $\mathcal{E}{\eta }^4$ is finite and in Eq. (3) we replace the $\sigma$ by its approximate value. Then it changes the estimation of error in calculation of the integral in order of $O\left(\frac{1}{N}\right)$.

In practice, when we simulate random variables $\xi$_, we receive simulation with some density $q\left(x\right)$ instead of simulation with the origin density $p\left(x\right)$. Now, we investigate the stability of the theoretical estimation Eq. (2). Let us consider the following expression:

$$\begin{array}{l}{\int }_{D}p\left(x\right)\frac{f\left(x\right)}{p\left(x\right)}\mathit{dx}-{\int }_{D}p\left(x\right)\frac{q\left(x\right)}{p\left(x\right)}\mathit{dx}={\int }_D\frac{f\left(x\right)}{p\left(x\right)}\left(p\left(x\right)-q\left(x\right)\right)\mathit{dx}\le \\ \le \epsilon {\int }_D\frac{\mid f\left(x\right)\mid }{p\left(x\right)}\mathit{dx},\end{array}$$E6

where $\epsilon =\underset{x\in D}{sup}\mid p\left(x\right)-q\left(x\right)\mid$. By $I\left(|f|/p\right)$ denote the integral $\int \frac{\mid f\left(x\right)\mid }{p\left(x\right)}\mathit{dx}$. The both values $\epsilon$ and $I\left(|f|/p\right)$ provide the guaranteed proximity of the real estimation to the theoretical one of the integral.

Example 1. The inequality Eq. (6) is reduced to the equality if $D=D_1\cup D_2,D_1\cap D_2=\varnothing$. For $x\in D_1$ we get $p\left(x\right)-q\left(x\right)\equiv \epsilon >0$, and $f\left(x\right)\equiv 0$ in $D_2$. If the condition $I\left(|f|/p\right)=+\infty$ holds, then the error of the real estimation of the integral will be infinity for any $\epsilon >0$.

Hence, a quality of pseudorandom variables (i.e., smallness of $\epsilon$) does not yet guaranties the smallness of the error in general, as the integral $I\left(|f|/p\right)$ have to be both finite and not very great in magnitude.

Suppose we simultaneously make the estimations for both $I\left(|f|/p\right)$ and the origin integral Eq. (2) using the same density $p\left(x\right)$. Then we need to ask boundedness of the integral ${\int }_D\frac{\mid f\left(x\right)\mid }{p^2\left(x\right)}\mathit{dx}$ to get the guaranteed stability of the estimation for the integral $I\left(|f|/p\right)$, and so on. The qualitative comparison of simulation with both densities $p_1\left(x\right)$ and $p_2\left(x\right)$ can be provided not only by a magnitude of the variance estimation (here we do not pay attention to the complexity of random values simulation) but also magnitudes of both the integrals $I\left(|f|/p_1\right)$ and $I\left(|f|/p_2\right)$. On the other hand we have the Schwarz inequality:

$${\int }_D\frac{f\left(x\right)}{p\left(x\right)}\left(p\left(x\right)-q\left(x\right)\right)\mathit{dx}\le \sqrt{{\int }_D{\left[p\left(x\right)-q\left(x\right)\right]}^2\mathit{dx}}\sqrt{{\int }_D\frac{f^2\left(x\right)}{p^2\left(x\right)}\mathit{dx}}.$$E7

The sufficient condition for the estimation to be stability is that ${\int }_D\frac{f^2\left(x\right)}{p^2\left(x\right)}$ to be finite and not very great. The Schwarz inequality is reduced to the equality if and only if

$$\lambda \frac{f\left(x\right)}{p\left(x\right)}=p\left(x\right)-q\left(x\right),$$E8

where $\lambda$ is a real number. From the above we get

$$\begin{array}{c}q\left(x\right)=p\left(x\right)-\lambda \frac{f\left(x\right)}{p\left(x\right)},\\ {\int }_{D}q\left(x\right)\mathit{dx}={\int }_{D}p\left(x\right)\mathit{dx}-\lambda {\int }_D\frac{f\left(x\right)}{p\left(x\right)}\mathit{dx}.\end{array}$$E9

Therefore, the equality in Eq. (7) is reached under the necessary condition ${\int }_D\frac{f\left(x\right)}{p\left(x\right)}\mathit{dx}=0$, when ${\int }_{D}p\left(x\right)\mathit{dx}={\int }_{D}q\left(x\right)\mathit{dx}=1$.

Example 2. The condition above is realized, e.g., if

$$\begin{array}{r}D=\left[-1,1\right],D_1=\left[-1,0\right],D_1=\left[0,1\right],\\ f\left(x\right)=-1\text{in}{D}_1,f\left(x\right)=1\text{in}{D}_2,p\left(x\right)=p\left(-x\right).\end{array}$$

Let us ${\eta }_k$ be ${\left[f\left(\xi \right)/p\left(\xi \right)\right]}^k$. Now we consider an estimation:

$$\mathcal{E}{\eta }_k={\int }_{D}p\left(x\right){\left[\frac{f\left(x\right)}{p\left(x\right)}\right]}^k\mathit{dx}={\int }_D\frac{f^k\left(x\right)}{p^{k-1}\left(x\right)}\mathit{dx},k\ge 1.$$E10

This expectation is actually estimated by the integral: ${\int }_{D}q\left(x\right){\left[f\left(x\right)/p\left(x\right)\right]}^k\mathit{dx}$:

$$\begin{array}{l}{\int }_D\frac{f^k\left(x\right)}{p^{k-1}\left(x\right)}\mathit{dx}-{\int }_{D}q\left(x\right)\frac{f^k\left(x\right)}{p^k\left(x\right)}\mathit{dx}={\int }_D\left[p\left(x\right)-q\left(x\right)\right]\frac{f^k\left(x\right)}{p^k\left(x\right)}\mathit{dx}\le \\ \le \epsilon {\int }_D\frac{\mid f^k\left(x\right)\mid }{p^k\left(x\right)}\mathit{dx}.\end{array}$$E11

The last integral is assumed to be a finite, and not very large. These conditions are desirable. In Eq. (11) the equality is reached like to the Example 1.

For all cases above, the stability will be observed if $\mid f\left(x\right)/p\left(x\right)\mid \le M<+\infty$ for not very great $M$. From the Schwarz inequality we have:

$$\begin{array}{l}{\int }_D\left[p\left(x\right)-q\left(x\right)\right]\frac{f^k\left(x\right)}{p^k\left(x\right)}\mathit{dx}={\int }_D\frac{p\left(x\right)-q\left(x\right)}{p^{\beta }\left(x\right)}\cdot \frac{f^k\left(x\right)}{p^{k-\beta }\left(x\right)}\mathit{dx}\le \\ \le \sqrt{{\int }_D\frac{{\left[p\left(x\right)-q\left(x\right)\right]}^2}{p^{2\beta }\left(x\right)}\mathit{dx}}\sqrt{{\int }_D\frac{f^{2k}\left(x\right)}{p^{2\left(k-\beta \right)}\left(x\right)}\mathit{dx}},\end{array}$$E12

where $\beta$ is a real number. We get a family of proximity measures for the distribution densities:

$${\int }_D\frac{{\left[p\left(x\right)-q\left(x\right)\right]}^2}{p^{2\beta }\left(x\right)}\mathit{dx}.$$E13

For $\beta =0,5$ we obtain expression

$${\chi }^2\left(p,q\right)={\int }_D\frac{{\left[p\left(x\right)-q\left(x\right)\right]}^2}{p^{(}x)}\mathit{dx}.$$E14

that is well known in the mathematical statistics.

For $\beta =-0,5$ we get ${\int }_D{\left[p\left(x\right)-q\left(x\right)\right]}^{2}p\left(x\right)\mathit{dx}$ and have the obvious inequalities

$$\begin{array}{l}{\int }_D{\left[p\left(x\right)-q\left(x\right)\right]}^{2}p\left(x\right)\mathit{dx}\le \underset{x\in D}{sup}p\left(x\right)\cdot {\int }_D{\left[p\left(x\right)-q\left(x\right)\right]}^2\mathit{dx},\\ {\int }_D{\left[p\left(x\right)-q\left(x\right)\right]}^{2}p\left(x\right)\mathit{dx}\le \underset{x\in D}{sup}{\left[p\left(x\right)-q\left(x\right)\right]}^2{\int }_{D}p\left(x\right)\mathit{dx}=\underset{x\in D}{sup}{\left[p\left(x\right)-q\left(x\right)\right]}^2.\end{array}$$

In Eq. (6) the equality is satisfied if and only if

$$\frac{p-q}{p^{\beta }}=\lambda \frac{f^k}{p^{k-\beta }},p-q=\lambda \frac{f^k{p}^{\beta }}{p^{k-\beta }},$$E15

i.e., the necessary condition is ${\int }_D\frac{f^k\left(x\right)}{p^{k-2\beta }\left(x\right)}\mathit{dx}=0.$ This is realized in the Example 2.

Let us remark that for $k=1$ and $\beta =1$ we have

$$\begin{array}{l}{\int }_D\frac{p\left(x\right)-q\left(x\right)}{p\left(x\right)}f\left(x\right)\mathit{dx}={\int }_D\frac{p-q}{\sqrt{p}}\cdot \frac{f}{\sqrt{p}}\le \\ \le \sqrt{{\int }_D\frac{{\left(p-q\right)}^2}{p}\mathit{dx}}\sqrt{{\int }_D\frac{f^2}{p}\mathit{dx}}=\sqrt{{\chi }^2\left(p,q\right)}\cdot \mathcal{E}{\eta }_2.\end{array}$$E16

We assume that integrals ${\int }_D\left[f\left(x\right)+{\epsilon }_i\left(x\right)\right]\mathit{dx}$ are known and the subintegral functions $f\left(x\right)+{\epsilon }_i\left(x\right)>0$ close to a function $f\left(x\right)\ge 0$. Suppose also

$$\begin{array}{l}1-{\delta }_1\left({\epsilon }_i\right)<\frac{f\left(x\right)}{f\left(x\right)+{\epsilon }_i\left(x\right)}<1+{\delta }_2\left({\epsilon }_i\right),\\ J\left(f\right)-{\delta }_3\left({\epsilon }_i\right)\le J\left(f+{\epsilon }_i\right)\le J\left(f\right)+{\delta }_4\left({\epsilon }_i\right),{\delta }_j\to 0,\mathrm{as}{\epsilon }_i\to 0,j=\mathrm{1,2,3,4}.\\ {\int }_{D}q\left(x\right)\frac{f\left(x\right)}{p\left(x\right)}\mathit{dx}={\int }_{D}q\left(x\right)\frac{f\left(x\right)J\left(f+{\epsilon }_i\right)}{f\left(x\right)+{\epsilon }_i\left(x\right)}\mathit{dx}=\mathcal{E}\widehat{\eta },\end{array}$$E17

where the random value $\widehat{\eta }$ has a form:

$$\widehat{\eta }=\frac{f\left(\widehat{\xi }\right)J\left(f+{\epsilon }_i\right)}{f\left(\widehat{\xi }\right)+{\epsilon }_i\left(\widehat{\xi }\right)}.$$E18

Here $\widehat{\xi }$ is distributed with the density $q\left(x\right)$. Keeping the above factors in mind we get the following:

$$\left[J\left(f\right)-{\delta }_3\right]\left[1-{\delta }_1\right]\le \eta \le \left[1+{\delta }_2\right]\left[J\left(f\right)+{\delta }_4\right];$$

1. Regardless of $q\left(x\right)$, i.e., regardless of quality of a pseudorandom number generator we have $\mathcal{E}\widehat{\eta }\to J\left(f\right),\mathit{var}\widehat{\eta }\to 0,$ as$\epsilon \to 0$;

2. All moments $\mathcal{E}{\widehat{\eta }}^k$ of the random variable are finite.

Let we calculate ${\int }_D\frac{\mid f\left(x\right)\mid }{p\left(x\right)}\mathit{dx}$ using the density

$$p_1\left(x\right)=\frac{\mid f\left(x\right)\mid }{I\left(|f|/p\right)p\left(x\right)},$$

then the estimation variance equals to zero.

Suppose we calculate the integral ${\int }_{D}f\left(x\right)\mathit{dx}$ with the density $p_1\left(x\right)$. In this case it would be interesting to know both the values ${\int }_D\frac{\mid f\left(x\right)\mid }{p_1}\mathit{dx}$ and ${\int }_D\frac{f^2}{p_1^2}\mathit{dx}$.

Proposition 1.$I\left(|f|/p_1\right)=I\left(|f|/p\right)$.

Proposition 2.$I\left(f^2/p_1^2\right)=I^2\left(|f|/p\right)I\left(p^2\right)$.

Now we consider the density

$$p_2\left(x\right)=\frac{f^2\left(x\right)}{p\left(x\right)I\left(f^2/p\right)}.$$

Using the density above for the estimation of the integral ${\int }_D\frac{f^2}{p}\mathit{dx}$ we obtain the estimation variance equals to zero.

Further we estimate the integral ${\int }_{D}f\left(x\right)\mathit{dx}$ with the density $p_2$: $\widehat{\eta }=f\left({\xi }_2\right)/p\left({\xi }_2\right)$, ${\xi }_2$ is distributed with $p_2\left(x\right)$. Suppose $\eta =f\left(\xi \right)/p\left(\xi \right)$, $\xi$ is distributed with $p\left(x\right)$; then $I\left(f^2/p_2\right)=I\left(f^2/p\right)$.

Proposition 3.$\mathit{var}\eta =\mathit{var}\widehat{\eta }$.

Proposition 4.${\left[\underset{D}{\int }\frac{f\left(x\right)}{p\left(x\right)}\mathit{dx}\right]}^2\le \left(\mathit{vol}D\right)\cdot \underset{D}{\int }\frac{f^2\left(x\right)}{p^2\left(x\right)}\mathit{dx}$, where $\mathit{vol}D$ is the volume of the domain $D$.

Proposition 5. If $\mathit{vol}D=1$, $p_3\left(x\right)=\frac{f^2\left(x\right)}{I\left(f^2\right)},{\eta }_3=\frac{f\left({\xi }_3\right)}{p\left({\xi }_3\right)},$ where ${\xi }_3$ is a random variable distributed with the density of $p_3\left(x\right)$; ${\eta }_4=\frac{f\left({\xi }_4\right)}{p\left({\xi }_4\right)},$ where ${\xi }_4$ is a random variable distributed with the density of $p\left(x\right)\equiv 1,$ then $\mathit{var}{\eta }_3=\mathit{var}{\eta }_4$.

In actual practice normalization constants are usually unknown for both $p_1\left(x\right)$ and $p_2\left(x\right)$. But using densities close to them we can get the approximate equalities in the Prepositions 1, 2, 3.

Now we consider

$$I\left(f\right)=\underset{{\left[0,1\right]}^{10}}{\int }{x}_1{x}_2\dots x_{10}{\mathit{dx}}_1\dots {\mathit{dx}}_{10},$$E19

where the integration domain $D={\left[0,1\right]}^{10}$ is the 10-dimensional unit cube, the subintegral function $f\left(x\right)$ is equals to $x_1{x}_2\dots x_{10}.$ To realize algorithms of the statistical modeling at a computer it is necessary to set a number $N$ of realizations for random variable $\eta =f\left(\xi \right)/p\left(\xi \right)$, where $\xi$ is distributed with the density $p\left(\xi \right)$. In fact we realize the discrete set of numbers ${\xi }_i,i=1,\dots ,N,$ which we can consider to be realizations of some distribution $q_N\left(x\right)$.

In all our numerical computations we use the pseudorandom number generator: generator type mt19937 [3]. For $s=10,p\left(x\right)\equiv 1$ we have

$$\begin{array}{c}I\left(f\right)={\left(1/2\right)}^{10}\approx 9,7656\cdot {10}^{-4},I\left(f/p\right)=I\left(f\right),\\ I\left(f^2/p\right)=I\left(f^2/p^2\right)={\left(1/3\right)}^{10}\approx 1,6935\cdot {10}^{-5},\\ \sigma =\sqrt{{\left(1/3\right)}^{10}-{\left(1/2\right)}^{20}}\approx \mathrm{3,998}\cdot {10}^{-3}.\end{array}$$

Table 1 shows the empirical estimations $\widehat{I}$ and $\widehat{\sigma }$ for $I$ and $\sigma$_, respectively. Taking $p\left(x_i\right)=3x_i^2$ over each coordinate we get

$N$	$\widehat{I}$	$\widehat{\sigma }$
1,000,000	$\mathrm{9,785}\cdot {10}^{-4}$	$\mathrm{3,991}\cdot {10}^{-3}$
9,000,000	$\mathrm{9,768}\cdot {10}^{-4}$	$\mathrm{4,001}\cdot {10}^{-3}$
81,000,000	$\mathrm{9,766}\cdot {10}^{-4}$	$\mathrm{3,997}\cdot {10}^{-3}$
100,000,000	$\mathrm{9,765}\cdot {10}^{-4}$	$\mathrm{3,995}\cdot {10}^{-3}$

Table 1.

The results of numerical calculations for the integral Eq. (19) with the uniform density.

$$\begin{array}{c}I\left(f^2/p\right)={\left(1/3\right)}^{10},I\left(f/p\right)=\infty ,I\left(f^2/p^2\right)=\infty ,\\ \sigma =\sqrt{{\left(1/3\right)}^{10}-{\left(1/2\right)}^{20}}.\end{array}$$

The value $\sigma$ is the same as one for $p\left(x\right)\equiv 1$.

For $N=100000000$ the computer code outputs an error because of machine zero divide. The reason of this event is $\eta ={\prod }_1^{10}(1/\left(3{\sqrt[3]{\alpha }}_i\right)$, where ${\alpha }_i$ are the pseudorandom numbers with the uniform density in $\left(0,1\right)$ [3]. If the formulas of random numbers simulation generate division by very small numbers then such formulas are one more source of the algorithms instability in the statistical modeling. As seen in Table 2 the value of integral $I\left(f\right)$ is successfully estimated, but the empirical estimations of $\widehat{\sigma }$ are sufficiently different from the theoretical value $\sigma$. This result is explained by the following: $\mathcal{E}{\eta }^4=\infty ,I\left(f/p\right)=\infty ,I\left(f^2/p^2\right)=\infty .$ Taking $p\left(x_i\right)=2,6x_i^{1,6}$ we obtain $\sigma \approx \mathrm{1,223}\cdot {10}^{-3}$, $I\left(f/p\right)\approx 0,67556$, the finite value of $\mathcal{E}{\eta }^4$, and $I\left(f^2/p^2\right)=\infty$. Although the last estimation is infinite, but Table 3 shows that both values $I\left(f\right)$ and $\widehat{\sigma }$ are successfully calculated. The value of $\widehat{\sigma }$ is very close to $\sigma$.

$N$	$\widehat{I}$	$\widehat{\sigma }$
1,000,000	$\mathrm{9,766}\cdot {10}^{-4}$	$\mathrm{3,276}\cdot {10}^{-3}$
9,000,000	$\mathrm{9,754}\cdot {10}^{-4}$	$\mathrm{3,619}\cdot {10}^{-3}$
81,000,000	$\mathrm{9,765}\cdot {10}^{-4}$	$\mathrm{3,759}\cdot {10}^{-3}$
100,000,000	Inf	nan

Table 2.

The results of numerical calculations for the integral Eq. (19) with the density $p\left(x_i\right)=3x_i^2$.

$N$	$\widehat{I}$	$\widehat{\sigma }$
1,000,000	$\mathrm{9,769}\cdot {10}^{-4}$	$\mathrm{1,212}\cdot {10}^{-3}$
9,000,000	$\mathrm{9,760}\cdot {10}^{-4}$	$\mathrm{1,220}\cdot {10}^{-3}$
81,000,000	$\mathrm{9,766}\cdot {10}^{-4}$	$\mathrm{1,223}\cdot {10}^{-3}$

Table 3.

The results of numerical calculations for the integral Eq. (19) with the density $p\left(x_i\right)=2,6x_i^{1,6}.$

3. Integrals of very large dimensions

We are coming now to the question of calculation of an integral

$$I\left(f\right)=\underset{0}{\overset{\infty }{\int }}\dots \underset{0}{\overset{\infty }{\int }}{e}^{-\left(x_1+x_2+\dots +x_s\right)}{\mathit{dx}}_1{\mathit{dx}}_2\dots {\mathit{dx}}_{10}=1$$E20

with the distribution density $p\left(x\right)={\lambda }^s{e}^{-\lambda \left(x_1+x_2+\dots +x_s\right)}$. For $\lambda \ge 2$_, the estimation variance $\eta$ is infinity. For $0<\lambda <2$_, the variance will be finite. For $\lambda >1$_, we obtain $I\left(f/p\right)=\infty$ and $I\left(f^2/p^2\right)=\infty$. However, as seen in Table 4, the results of calculations for $N=10000$ allow us to make the conclusion below. If we have the pseudorandom generator of the high quality and a good $p\left(x\right)$ then we can calculate the very high dimensional integrals.

$\lambda$	$s$	$\widehat{I}$	$\widehat{\sigma }$	$\sigma$
1,01	1000	0,999	0,320	0,324
1,01	10,000	0,995	1,39	1,31
1005	20,000	1004	0,807	0,805
1003	40,000	0,994	0,649	0,658
1001	80,000	0,991	0,605	0,661

Table 4.

The results of numerical calculations for the integral Eq. (20).

4. Special integrals

Let us consider the following class of the integrals:

$$\underset{{\left[0,+\infty \right)}^s}{\int }f\left(x_1,x_2,\dots ,x_s\right){\mathit{dx}}_1{\mathit{dx}}_2\dots {\mathit{dx}}_s,$$E21

We define the behavior of the subintegral function as follows: $f\left(x_1,x_2,\dots ,x_s\right)$ to be [label = ()]

1. close to 1 in the cube ${\left[0,a\right]}^s$ as $0<a<1$;

2. much less than unit as $a<x_i<b,b\ge a$;

3. equil to zero as $b\le x_i<\infty$.

Note that very often the integration of functions can be reduced to the linear combination of the integrals similar to Eq. (21) using various replacements of variables.

Below, let us perform a theoretical and numerical analysis how to integrate a model function from our class. The model function is assumed to be $f\equiv 1$ as $0\le x_i\le a$, otherwise $f=0$. We take both distribution densities set $p_1\left(x_i\right)=\lambda e^{-\lambda x_i},0\le x_i<\infty$ and $p_2\left(x_i\right)=\left(\omega +1\right){\left(1-x_i\right)}^{\omega },0\le x_i<1$ to be examined. Our goal is to determine what of two densities provides the best accuracy of the integral computation with given model function.

If we simulate a random point $\xi =\left({\xi }_1,\dots ,{\xi }_s\right)$ with densities $p_1\left(x_i\right)=\lambda e^{-\lambda x_i}$ then the integral estimation is given by

$$\begin{array}{c}{\eta }_s=\prod _{i=1}^s\lambda e^{{\mathit{\lambda \xi }}_i},\\ \mathbf{var}{\eta }_s=\mathcal{E}{\eta }_s^2-{\left(\mathcal{E}{\eta }_s\right)}^2={\left[\frac{1}{\lambda }\underset{0}{\overset{a}{\int }}{e}^{\mathit{\lambda x}}\mathit{dx}\right]}^s-a^s={\left[\frac{e^{\mathit{\lambda a}}-1}{{\lambda }^2}\right]}^s-a^s.\end{array}$$E22

Testing the variance $\mathbf{var}{\eta }_s$ for the extremum over $\lambda$ we get the minimum condition

$$\lambda {\mathit{ae}}^{\mathit{\lambda a}}-2e^{\mathit{\lambda a}}+2=0.$$E23

Let $A$ be $\mathit{\lambda a}$; then the equation Eq. (23) is reduced to

$${\mathit{Ae}}^A-2e^A+2=0.$$E24

The equation above has the unique root at $A\approx 1,593620$. It follows that ${\lambda }_{\mathit{min}}=A/a$. For such ${\lambda }_{\mathit{min}}$ the relative error with the 3$\sigma$ rule is given by

$$\frac{3\sigma }{\sqrt{N}{a}^s}=3{\left(\frac{{\left(e^A-1\right)}^s}{A^{2s}}-1\right)}^{0,5}/\sqrt{N}.$$E25

Suppose $N=9\cdot {10}^6,s=10,\lambda =A/a$; then the theoretical value of the relative error is approximately $8,72\cdot {10}^{-3}$. The numerical estimation of the relative error is approximately $8,69\cdot {10}^{-3}$ as $a\in \left[0,1;\mathrm{0,001}\right]$. Thus, the numerical estimation of one gives a good fit to the predicted value over a wide range of $a$.

Now we discuss the use of the density $p_2\left(x_i\right)=\left(\omega +1\right){\left(1-x_i\right)}^{\omega }$. First, we estimate the second moment of a random value ${\eta }_s$:

$$\mathcal{E}{\eta }_s^2={\left[\frac{1-{\left(1-a\right)}^{1-\omega }}{1-{\omega }^2}\right]}^s.$$E26

The parameter $\omega$ is chosen to be $A/a$; then the expression above is rewritten as follows

$$\mathcal{E}{\eta }_s^2={\left[\frac{1}{1-A^2/a^2}\left(1-{\left(1-a\right)}^{1-A/a}\right)\right]}^s.$$E27

Let us consider $\mathcal{E}{\eta }_s^2$ as $a\to 0$:

$$\begin{array}{l}\underset{a\to 0}{lim}\mathcal{E}{\eta }_s^2=\underset{a\to 0}{lim}{\left[\frac{a^2}{a^2-A^2}\left(1-{\left(1-a\right)}^{1-A/a}\right)\right]}^s=\\ =-\underset{a\to 0}{lim}{\left[\frac{a^2}{A^2}\left(1-\left(1-a\right){\left(1-a\right)}^{-A/a}\right)\right]}^s=-\underset{a\to 0}{lim}{\left[\frac{a^2}{A^2}\left(1-\left(1-a\right){\left[{\left(1-a\right)}^{-1/a}\right]}^A\right)\right]}^s=\\ =-\underset{a\to 0}{lim}{\left[\frac{a^2}{A^2}\left(1-\left(1-a\right)e^A\right)\right]}^s=-\underset{a\to 0}{lim}{\left[\frac{a^2}{A^2}\left(1-e^A\right)\right]}^s\sim {\left[\frac{\left(e^A-1\right)}{A^2}{a}^2\right]}^s.\end{array}$$E28

Comparison between Eqs. (22) and (28) allows to make the following conclusion. If $w$ is chosen to be $A/a$ then the asymptotics of variances, as $a\to 0$, are the same in the densities set of $p_1\left(x_i\right),p_2\left(x_i\right)$. In the numerical simulation the relative accuracy of $\approx 8,68\cdot {10}^{-3}$ is reached as $N=9\cdot {10}^6,s=10,\omega =A/a,a\in \left[0,01;\mathrm{0,001}\right]$.

Let us turn now to the integral

$$I\left(f\right)=\underset{{\left[0,1\right]}^{10}}{\int }f\left(x_1,x_2,\dots ,x_{10}\right){\mathit{dx}}_1{\mathit{dx}}_2\dots {\mathit{dx}}_{10},$$E29

where

$$f\left(x_1,x_2,\dots ,x_{10}\right)=\left\{\begin{array}{ll}1,& 0\le x_i\le 1/4,\\ 0,& \text{otherwise}.\end{array}\right.$$E30

We put $p\left(x\right)\equiv 1$; then $I\left(f\right)=I\left(f/p\right)=I\left(f^2/p\right)={\left(1/4\right)}^{10}\approx 9,5367\cdot {10}^{-7}.$ For $N=640000$ all realizations are turned out to be equal to zero, i.e., $q_{640000}\left(x\right)=0$ as $0\le x_i\le 1/4$. In this case, we have

$$\underset{{\left[0,1\right]}^{10}}{\int }p\left(x\right)\frac{f\left(x\right)}{p\left(x\right)}\mathit{dx}-\underset{{\left[0,1\right]}^{10}}{\int }{q}_{640000}\frac{f\left(x\right)}{p\left(x\right)}\mathit{dx}={\left(1/4\right)}^{10}-0={\left(1/4\right)}^{10}.$$E31

In accordance with $N$_, the realizations numbers of ${\eta }_i$ are turned out to be equal to 1 as $N=810000$; equal to 3 as $N=4000000$; equal to 15 as $N=16000000$.

We now take $p\left(x_i\right)=\left(\omega +1\right){\left(1-x_i\right)}^{\omega }$ in the unit cube ${\left[0,1\right]}^{10}\left(\omega >1\right)$. Such choice provides the gross realizations of points in ${\left[0,1/4\right]}^{10}$ and as consequence, we get benefit in quality of random values (simultaneously, we have decrease of the estimation variance, and as consequence decrease of the statistical error with the 3$\sigma$ rule.) Table 5 shows the calculations results for $N=9000000$. Note that $\widehat{\sigma }$ reaches the minimum as $\omega =5$. In this case, we have $\widehat{p}\left(x\right)=6{\left(1-x\right)}^5,I\left(f/\widehat{p}\right)=I\left(f^2/\widehat{p}\right)\approx 3,49\cdot {10}^{-11}.$ Making the more detailed research for both $\omega =5$ and the theoretical value $\sigma \approx 5,83\cdot {10}^{-6}$ we get the results represented in Table 6. If the function $f\left(x_1,x_2,\dots ,x_{10}\right)$ close to some constant in ${\left[0,1/4\right]}^{10}$ and small out of this interval then we can advise to use $\widehat{p}=6{\left(1-x\right)}^5$ to calculate the integral in ${\left[0,1\right]}^{10}$.

$w$	$\widehat{I}$	rule “$3\sigma$”
3	$9,53\cdot {10}^{-7}$	$8,61\cdot {10}^{-9}$
4	$9,58\cdot {10}^{-7}$	$6,36\cdot {10}^{-9}$
5	$9,55\cdot {10}^{-7}$	$5,82\cdot {10}^{-9}$
6	$9,55\cdot {10}^{-7}$	$6,49\cdot {10}^{-9}$
7	$9,51\cdot {10}^{-7}$	$8,06\cdot {10}^{-9}$

Table 5.

The results of numerical calculations for the integral Eq. (29) at various $\omega$ values.

$N$	$\widehat{I}$	$\widehat{\sigma }$
640,000	$9,61\cdot {10}^{-7}$	$6,15\cdot {10}^{-6}$
810,000	$9,56\cdot {10}^{-7}$	$6,00\cdot {10}^{-6}$
1,000,000	$9,51\cdot {10}^{-7}$	$5,91\cdot {10}^{-6}$
4,000,000	$9,52\cdot {10}^{-7}$	$5,74\cdot {10}^{-6}$
16,000,000	$9,55\cdot {10}^{-7}$	$5,84\cdot {10}^{-6}$

Table 6.

The results of numerical calculations for the integral Eq. (29) at $\omega =5$.

5. Moments calculation

We are now concerned with the following issue: to find the kth moments of a random value $\tau$ with the distribution density $p\left(x\right)$:

$$\mathcal{E}{\tau }^k=\underset{a}{\overset{b}{\int }}{x}^{k}p\left(x\right)\mathit{dx}.$$E32

In fact we have realizations of the random value $\xi$ with a distribution density $q\left(x\right)$. With $p\left(x\right)$ replaced by $q\left(x\right)$ in Eq. (32) we get an error

$$\underset{a}{\overset{b}{\int }}{x}^{k}p\left(x\right)\mathit{dx}-\underset{a}{\overset{b}{\int }}{x}^{k}q\left(x\right)\mathit{dx}=\underset{a}{\overset{b}{\int }}{x}^k\left[p\left(x\right)-q\left(x\right)\right]\mathit{dx}.$$E33

Suppose $b=\infty$ and ${\xi }_{\mathit{max}}$ are the maximum value of the random variable over the all realizations for fixed $N$; then value of ${\xi }_{\mathit{max}}$ gives shift ${\int }_{{\xi }_{\mathit{max}}}^{\infty }{x}^{k}p\left(x\right)\mathit{dx}$ that increases both monotonically and without limit. The condition $q\left(x\right)=0$ as $x>{\xi }_{\mathit{max}}$ determines the lower limit of the last integral.

Many solutions of the boundary-value problems for the elliptic and parabolic Equations [4, 5] have a form of the expectations for the random value moments. Meaning of these expectations is the first exit time of the Wiener process trajectories to the domain boundary.

Let a domain be the three-dimensional ball with the radius $r=1$ and the Wiener trajectories start from the ball center; then a function of distribution of the first exit time for the Wiener trajectory is, in particular, given by [5].

$$F\left(t\right)=1+2\sum _{k=1}^{\infty }{\left(-1\right)}^{k}exp\left(-k^2{\pi }^{2}t/2\right),t\in \left[0,+\infty \right).$$E34

From the above, we obtain the distribution density:

$$p\left(t\right)=2\sum _{k=1}^{\infty }{\left(-1\right)}^{k+1}\mu k^{2}exp\left(-\mu k^{2}t\right),\mu ={\pi }^2/2.$$E35

Assuming $\tau$ is distributed with this density and calculating the expectation of the $k$th moment we get

$$\mathcal{E}{\tau }^k=\underset{0}{\overset{\infty }{\int }}{t}^{k}p\left(t\right)\mathit{dt}.$$E36

In Table 7, we put the calculations results for $N=1000000$. The $k$th moment expectation can be represented in a form.

Moment	Simulation	Theory
1	$\mathrm{3,304}\cdot {10}^{-1}$	$\mathrm{3,333}\cdot {10}^{-1}$
2	$\mathrm{1,553}\cdot {10}^{-1}$	$\mathrm{1,556}\cdot {10}^{-1}$
3	$\mathrm{9,848}\cdot {10}^{-2}$	$\mathrm{9,841}\cdot {10}^{-2}$
4	$\mathrm{8,076}\cdot {10}^{-2}$	$\mathrm{8,063}\cdot {10}^{-2}$
5	$\mathrm{8,291}\cdot {10}^{-2}$	$\mathrm{8,193}\cdot {10}^{-2}$
6	$\mathrm{9,843}\cdot {10}^{-2}$	$\mathrm{9,969}\cdot {10}^{-2}$
7	$\mathrm{1,319}\cdot {10}^{-1}$	$\mathrm{1,414}\cdot {10}^{-1}$
8	$\mathrm{2,070}\cdot {10}^{-1}$	$\mathrm{2,293}\cdot {10}^{-1}$
9	$\mathrm{4,518}\cdot {10}^{-1}$	$\mathrm{4,182}\cdot {10}^{-1}$
10	$\mathrm{7,286}\cdot {10}^{-1}$	$\mathrm{8,474}\cdot {10}^{-1}$
11	$\mathrm{9,021}\cdot {10}^2$	$\mathrm{4,251}\cdot {10}^3$
12	$\mathrm{1,183}\cdot {10}^3$	$\mathrm{1,637}\cdot {10}^4$
13	$\mathrm{8,389}\cdot {10}^3$	$\mathrm{6,634}\cdot {10}^4$

Table 7.

The results of numerical calculations for the moments by the first way.

$$\mathcal{E}{\tau }^k=\underset{0}{\overset{\infty }{\int }}q\left(x\right)t^k\frac{p\left(t\right)}{q\left(t\right)}\mathit{dt},$$E37

where $q\left(t\right)$ is some density in $\left[0,\infty \right)$. Taking $q\left(t\right)=\lambda exp\left(-\mathit{\lambda t}\right),$ for $\lambda ={\pi }^2/2$ we get that the number of realizations ${\xi }_i>1$ will be almost twice as small as in the case of the modeling with the original $p\left(t\right)$. In this situation we should obtain degradation of the estimation for the high moments. The calculations results with $q\left(t\right)$ for $N=1000000$ are represented in Table 8. However, in realizations at a computer we get the obvious improvement in quality of the moments estimation for all $k$ from 5 to 20. Consider the choice of modeling strategy with regards to the variance. Suppose $\xi$ and $\eta$ be estimations of the statistical modeling for a value $J$, i.e., $\mathcal{E}\xi =\mathcal{E}\eta =J$ with the variances of ${\sigma }_1^2\left(\xi \right),{\sigma }_2^2\left(\eta \right)$ and the realizations of $\xi =\left({\xi }_1+\dots +{\xi }_N\right)/N$, $\eta =\left({\eta }_1+\dots +{\eta }_N\right)/N$. It would seem that for ${\sigma }_1\left(\xi \right)<{\sigma }_2\left(\eta \right)$ the real estimation of $\xi$ will be occurred close to the origin value of $J$. But this statement does not need to be always true. Without loss of generality it can believed that $J=0$. Additionally, if $N$ is large enough then $\xi$, $\eta$ are chosen be normal random variables with $N\left(0,{\sigma }_1\right)$ and $N\left(0,{\sigma }_2\right)$, respectively. The following theorem holds.

Moment	Simulation
5	$\mathrm{8,188}\cdot {10}^{-2}$
6	$\mathrm{1,013}\cdot {10}^{-1}$
7	$\mathrm{1,468}\cdot {10}^{-1}$
8	$\mathrm{2,247}\cdot {10}^{-1}$
9	$\mathrm{3,950}\cdot {10}^{-1}$
10	$\mathrm{8,283}\cdot {10}^{-1}$
18	$\mathrm{2,833}\cdot {10}^3$
19	$\mathrm{1,056}\cdot {10}^4$
20	$\mathrm{5,118}\cdot {10}^4$

Table 8.

The results of numerical calculations for the moments by the second way.

Proposition 6. Let $\xi ,\eta$ be normal random variables, and $\xi \sim N\left(0,{\sigma }_1\right),\eta \sim N\left(0,{\sigma }_2\right)$ then $P\left(\xi |>|\eta |\right)=\frac{2}{\pi }arctan\frac{{\sigma }_1}{{\sigma }_2}.$

Proof:

$$\begin{array}{l}P|\xi >\left|\eta \right|=\frac{1}{{\sigma }_1\sqrt{2\pi }}\underset{-\infty }{\overset{0}{\int }}{e}^{-\frac{y^2}{2{\sigma }_1^2}}\mathit{dy}\left\{\frac{1}{{\sigma }_2\sqrt{2\pi }}\underset{0}{\overset{y}{\int }}{e}^{-\frac{x^2}{2{\sigma }_2^2}}\mathit{dx}+\frac{1}{{\sigma }_2\sqrt{2\pi }}\underset{0}{\overset{\left|y\right|}{\int }}{e}^{-\frac{x^2}{2{\sigma }_2^2}}\mathit{dx}\right\}+\\ +\frac{1}{{\sigma }_1\sqrt{2\pi }}\underset{0}{\overset{\infty }{\int }}{e}^{-\frac{y^2}{2{\sigma }_1^2}}\mathit{dy}\left\{\frac{1}{{\sigma }_2\sqrt{2\pi }}\underset{0}{\overset{y}{\int }}{e}^{-\frac{x^2}{2{\sigma }_2^2}}\mathit{dx}+\frac{1}{{\sigma }_2\sqrt{2\pi }}\underset{0}{\overset{\left|y\right|}{\int }}{e}^{-\frac{x^2}{2{\sigma }_2^2}}\mathit{dx}\right\}=\\ =\frac{2}{{\sigma }_1\sqrt{2\pi }}\underset{0}{\overset{\infty }{\int }}{e}^{-\frac{y^2}{2{\sigma }_1^2}}\mathit{dy}\left\{\frac{2}{{\sigma }_2\sqrt{2\pi }}\underset{0}{\overset{y}{\int }}{e}^{-\frac{x^2}{2{\sigma }_2^2}}\mathit{dx}\right\}=\frac{4}{2{\mathit{\pi \sigma }}_1{\sigma }_2}\underset{0}{\overset{\infty }{\int }}{e}^{-\frac{y^2}{2{\sigma }_1^2}}\mathit{dy}\cdot \underset{0}{\overset{y}{\int }}{e}^{-\frac{x^2}{2{\sigma }_2^2}}\mathit{dx}.\end{array}$$

Using Taylor expansion

$$e^{-\frac{y^2}{2{\sigma }_2^2}}=\sum _{n=0}^{\infty }{\left(-1\right)}^n\frac{y^{2n}}{2^n{\sigma }_2^{2n}n!},$$

we get

$$\begin{array}{l}\underset{0}{\overset{y}{\int }}\sum _{n=0}^{\infty }{\left(-1\right)}^n\frac{x^{2n}}{2^n{\sigma }_2^{2n}n!}\mathit{dx}=\sum _{n=0}^{\infty }{\left(-1\right)}^n\frac{y^{2n+1}}{\left(2n+1\right)2^n{\sigma }_2^{2n}n!},\\ \underset{0}{\overset{\infty }{\int }}\sum _{n=0}^{\infty }{\left(-1\right)}^n\frac{y^{2n+1}}{\left(2n+1\right)2^n{\sigma }_2^{2n}n!}{e}^{-\frac{y^2}{2{\sigma }_1^2}}\mathit{dy}=\sum _{n=0}^{\infty }{\left(-1\right)}^n\frac{1}{\left(2n+1\right)2^n{\sigma }_2^{2n}n!}\underset{0}{\overset{\infty }{\int }}{y}^{2n+1}{e}^{-\frac{y^2}{2{\sigma }_1^2}}\mathit{dy}=\\ =\sum _{n=0}^{\infty }{\left(-1\right)}^n\frac{1}{\left(2n+1\right)2^n{\sigma }_2^{2n}n!}\cdot \frac{2^{n+1}{\sigma }_1^{2n+2}n!}{2}.\left(\ast \right)\end{array}$$

Note that the last equality is obtained with the help of the formula:

$$\underset{0}{\overset{\infty }{\int }}{x}^{2n+1}{e}^{-{\mathit{px}}^2}\mathit{dx}=\frac{n!}{2p^{n+1}},p>0.$$

In our case, $p$ is $\frac{1}{2{\sigma }_1^2}$. We continue the equalities chain which is broken at (*):

$$\begin{array}{l}\left(\ast \right)=\sum _{n=0}^{\infty }{\left(-1\right)}^n\frac{1}{\left(2n+1\right)}\frac{{\sigma }_1^{2n+2}}{{\sigma }_2^{2n}}=\sum _{n=0}^{\infty }{\left(-1\right)}^n\frac{{\sigma }_1^2}{\left(2n+1\right)}{\left(\frac{{\sigma }_1}{{\sigma }_2}\right)}^{2n}=\\ =\frac{2}{\pi }\sum _{n=0}^{\infty }{\left(-1\right)}^n\frac{{\left({\sigma }_1/{\sigma }_2\right)}^{2n+1}}{2n+1}=\frac{2}{\pi }arctan\frac{{\sigma }_1}{{\sigma }_2}.\end{array}$$

For the k-moment calculation we take

$$q_k\left(t\right)=\frac{{\lambda }^{k+1}{t}^k{e}^{-\mathit{\lambda t}}}{k!},\lambda =\frac{{\pi }^2}{2}$$

and get the results shown in Table 9.

Moment	Simulation
1	$\mathrm{3,297}\cdot {10}^{-1}$
2	$\mathrm{1,556}\cdot {10}^{-1}$
3	$\mathrm{9,842}\cdot {10}^{-2}$
4	$\mathrm{8,066}\cdot {10}^{-2}$
5	$\mathrm{8,194}\cdot {10}^{-2}$
6	$\mathrm{9,968}\cdot {10}^{-2}$
7	$\mathrm{1,414}\cdot {10}^{-1}$
8	$\mathrm{2,293}\cdot {10}^{-1}$
9	$\mathrm{4,181}\cdot {10}^{-1}$
10	$\mathrm{8,474}\cdot {10}^{-1}$
18	$\mathrm{4,251}\cdot {10}^3$
19	$\mathrm{1,637}\cdot {10}^4$
20	$\mathrm{6,634}\cdot {10}^4$

Table 9.

The results of numerical calculations for the moments by the third way.

6. Integral from the regression function

Now, we consider the issue of calculation of an integral

$$\underset{D}{\int }f\left(x\right)\mathit{dx},$$E38

where the function $f\left(x\right)$ has no an analytical expression. Suppose there exists a random variable $\xi \left(x,w\right)$ such that its expectation is equals to $\mathcal{E}\xi \left(x,w\right)=f\left(x\right)$ for some fixed $x$. The random variable $\xi \left(x,w\right)$ may be realized neither as result of the physical measurements or some calculations (e.g., using the modeling statistical method). In this case the optimal density is given by [6].

$$p\left(x\right)=\frac{f\left(x\right)}{\sqrt{d\left(x\right)+\lambda }},$$E39

where $d\left(x\right)$ is the variance of the random variable $\xi \left(x,w\right)$. Note that one should use the optimal density from [1] if complexity in calculations (experimental measurements) is much different from each other for any $x$. We determine the parameter $\lambda$ from the condition $\underset{D}{\int }f\left(x\right)/\sqrt{d\left(x\right)+\lambda }\mathit{dx}=1.$ Really in practice, we find a priori or a posteriori approaches to both $f\left(x\right)$ and $d\left(x\right)$. By $\overline{f}\left(x\right)$ and $\overline{d}\left(x\right)$ denote these approaches. Then the approach to the optimal $p\left(x\right)$ will look like

$$\overline{p}\left(x\right)=\frac{\overline{f}\left(x\right)}{\sqrt{\overline{d}\left(x\right)+\overline{\lambda }}}\text{and}\underset{D}{\int }\frac{\overline{f}\left(x\right)}{\sqrt{\overline{d}\left(x\right)+\overline{\lambda }}}\mathit{dx}=1.$$E40

The parameter $\lambda$ is often turned out to be find enough complicity [6]. If the domain $D$ is the interval $\left[0,H\right]$ for small $H$ then it is suppose to use the quasioptimal density $\overline{p}\left(x\right)$.

Example 3. We now consider the following issue: Suppose $f\left(x\right)=x$, $d\left(x\right)=1/x,D=\left[0,H\right]$. The optimal density is given by

$$p\left(x\right)=\frac{x}{\sqrt{1/x+\lambda }}=\frac{x\sqrt{x}}{\sqrt{\lambda \sqrt{x}+1}}\sim c\cdot x^{3/2}.$$E41

We take the quasioptimal density in the form $\overline{p}\left(x\right)=5H^{5/2}{x}^{3/2}/2.$ In this case, for $p\left(x\right)\equiv 1$ the estimation variance of the random value $\eta =f\left(x,w\right)/p\left(x\right)$:

$$\mathbf{var}\eta =\underset{0}{\overset{H}{\int }}d\left(x\right)p\left(x\right)\mathit{dx}+\underset{0}{\overset{1}{\int }}\frac{f^2}{p\left(x\right)}\mathit{dx}-I^2$$E42

is equals to $\infty$. Taking $p\left(x\right)=2x/H^2$ we get $\mathbf{var}\eta =2/H.$ But if the function $f\left(x\right)$ was precisely known for the same density $p\left(x\right)$ then $\mathbf{var}\eta =0.$ If we choose the quasioptimal density $\overline{p}\left(x\right)=5H^{5/2}{x}^{3/2}/2$ then the estimation variance of $\eta$ is equals to $17H^4/12+4/\left(15H\right)$. For $H\to 0$ the variance behaves approximately as $4/\left(15H\right)$. It is much the better than $2/H$. For $H=1$ the estimation variance with the density $p\left(x\right)=2x/H^2$ is equals to 2, and the estimation variance with the quasioptimal density $\overline{p}\left(x\right)$ is equals to 101/60.

Suppose we practically realize calculation of the integral Eq. (38) with $d\left(x\right)=1/x$; then one should discard the interval $\left[0,\delta \right]$ and to calculate ${\int }_{\delta }^{H}f\left(x\right)\mathit{dx}$ because of the values $\mid \xi \left(x,w\right)\mid$ can be the intolerably large. Also one should replace $f\left(x\right)$ by $\widehat{f}\left(x\right)$:

$$\widehat{f}\left(x\right)=\left\{\begin{array}{ll}0,& 0\le x\le \delta ,\\ x,& \delta <x\le H.\end{array}\right.$$E43

The shift is $\underset{0}{\overset{\delta }{\int }}x\mathit{dx}={\delta }^2/2$ and choosing $\delta \sim 1/\sqrt[4]{N}$ we get the total error ${\delta }^2/2+3\sigma /\sqrt{N}$ of oder $O\left(1/\sqrt{N}\right)$.

In applications the estimation variance for the integral functionals (e.g., field flow calculation neither across the arc or the surface) from the solutions of the boundary-value problems for both the linear [1] or nonlinear [2] elliptic equations is of interest. For the above variance is $d\left(x\right)\sim B/x^2,f\left(x\right)\approx a_0+a_{1}x+a_2{x}^2+\dots$, where $x$ is the distance to the domain boundary. Suppose $f\left(x\right)\approx a_{1}x+a_2{x}^2+\dots$; then the optimal density is given by

$$p\left(x\right)=\frac{a_{1}x+a_2{x}^2+\dots }{\sqrt{B/x^2+\lambda }}.$$E44

The quasioptimal density has the form $\overline{p}\left(x\right)=3x^2/H^3$ for small $H$ in $\left[0,H\right]$. In applications, this case is of our main interest. Taking $\delta =1/\sqrt[4]{N}$ like in the Example 3 we get the asymptotics of decrease for the total error as $O\left(1/\sqrt{N}\right)$.

Suppose $d\left(x\right)\sim B/x^2,f\left(x\right)\approx a_0+a_{1}x+a_2{x}^2+\dots$, and $a_0\ne 0$ then there is no density kind of $\overline{p}\left(x\right)=\left(w+1\right)x^w,x\in \left[0,H\right]$ with the finite variance. The density $p\left(x\right)=\mid f\left(x\right)\mid /\sqrt{d\left(x\right)+\lambda }$ will be give the estimation with the infinity variance. Instead of calculation of the integral $\underset{0}{\overset{H}{\int }}f\left(x\right)\mathit{dx}$ we will be calculate the integral $\underset{\delta }{\overset{H}{\int }}f\left(x\right)\mathit{dx}$. For this integral we already can choice the quasioptimal density with the finite variance of the estimation: $\overline{p}\left(x\right)=2x/\left(H^2-{\delta }^2\right)$. For $\delta \sim O\left(lnN/\sqrt{N}\right)$ the total error will have the asymptotics $O\left(lnN/\sqrt{N}\right)$.

Example 4. Suppose $d\left(x\right)=1/x^2,f\left(x\right)=1,H=1;$ then the asymptotics of the variance with the quasioptimal density has kind of $\left(-2,5\cdot ln\delta \right)$.

In conditions of Example 4, choice of the optimal density in the form

$$p\left(x\right)=\frac{x}{\sqrt{1-{\delta }^2}\sqrt{1-x^2}},x\in \left[\delta ,1\right]$$E45

yields the following result: the estimation variance will have asymptotics $\left(-2\cdot ln\delta \right)$ for $\delta \to 0$.

Remark. If we know that a value of $f\left(x\right)$ in the interval $\left[0,\delta \right]$ close to the number $f_0$, then in Eq. (43) to use

$$\widehat{f}\left(x\right)=\left\{\begin{array}{ll}{f}_0,& 0\le x\le \delta ,\\ x,& \delta <x\le H,\end{array}\right.$$E46

more efficiently and also to take $\underset{0}{\overset{H}{\int }}f\left(x\right)\mathit{dx}\approx f_0\delta +\underset{\delta }{\overset{H}{\int }}f\left(x\right)\mathit{dx}.$

7. Conclusion

In the paper we describe the sufficient conditions of the stable calculations for the multidimensional integrals by the Monte Carlo method. We get the results of numerous numerical computations using the mt19937 pseudorandom number generator. The article results can be also useful in the practical solution of the boundary value problem, for both the elliptic and parabolic equations. The earlier suggested approach to the optimal choice of the density [1, 6] often needs to solve a complicated secondary task. In the paper we suggest the approach to choice of the quasioptimal densities that is of considerable interest in applied problems solution.