**2. Sufficient conditions of stability**

Let

$$I = \int\_{D} f(\mathbf{x})d\mathbf{x} \tag{1}$$

be the Riemann integral. Here *D* is a domain of the s-dimensional Euclidean space *Rs* . 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 *η* ¼ *f*ð Þ*ξ =p*ð Þ*ξ* :

$$\int\_{D} f(\mathbf{x})d\mathbf{x} = \int\_{D} p(\mathbf{x}) \frac{f(\mathbf{x})}{p(\mathbf{x})} d\mathbf{x} = \mathcal{E} \frac{f(\mathbf{x})}{p(\mathbf{x})} = \mathcal{E}\eta. \tag{2}$$

Here *p x*ð Þ is a density of random variable *ξ*∈ *D*. We put *p x*ð Þ 6¼ 0 for *f x*ð Þ 6¼ 0, and we say that there exists integral

$$\int\_{D} |f(\mathbf{x})|d\mathbf{x}.$$

A variance of the random value *η*:

$$\sigma^2 = \mathbf{var}\eta = \mathcal{E}\eta^2 - \left(\mathcal{E}\eta\right)^2 = \int\_D p(\mathbf{x}) \left[\frac{f(\mathbf{x})}{p(\mathbf{x})}\right]^2 d\mathbf{x} - I^2 = \int\_D \frac{f^2(\mathbf{x})}{p(\mathbf{x})} d\mathbf{x} - I^2. \tag{3}$$

We estimate the mathematical expectation E*η* by the sum P*<sup>N</sup> <sup>i</sup>*¼<sup>1</sup>*ηi=N*, where *<sup>η</sup><sup>i</sup>* are independent realizations of the random value *η*. Suppose *σ*<sup>2</sup> is finite, and *N* is large enough; then from the classical central limit theorem, it follows that the random value P*<sup>N</sup> <sup>i</sup>*¼<sup>1</sup>*ηi=<sup>N</sup>* has distribution close to the normal distribution with a mathematical expectation *I*, and mean-square deviation *σ=* ffiffiffiffi *N* <sup>p</sup> . This property above is useful to estimate error, e.g., using the 3*σ* rule. So we have

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

with probability 0,997, approximately.

Let us *η<sup>i</sup>* be realizable sampling values; then the value *σ* 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. \tag{5}
$$

Suppose E*η*<sup>4</sup> is finite and in Eq. (3) we replace the *σ* by its approximate value. Then it changes the estimation of error in calculation of the integral in order of *O* <sup>1</sup> *N* � �.

In practice, when we simulate random variables *ξ*, we receive simulation with some density *q x*ð Þ instead of simulation with the origin density *p x*ð Þ. Now, we investigate the stability of the theoretical estimation Eq. (2). Let us consider the following expression:

$$\begin{split} \int\_{D} p(\mathbf{x}) \frac{f(\mathbf{x})}{p(\mathbf{x})} d\mathbf{x} - \int\_{D} p(\mathbf{x}) \frac{q(\mathbf{x})}{p(\mathbf{x})} d\mathbf{x} &= \int\_{D} \frac{f(\mathbf{x})}{p(\mathbf{x})} (p(\mathbf{x}) - q(\mathbf{x})) d\mathbf{x} \le \\ &\le e \int\_{D} \frac{|f(\mathbf{x})|}{p(\mathbf{x})} d\mathbf{x}, \end{split} \tag{6}$$

*Stability of Algorithms in Statistical Modeling DOI: http://dx.doi.org/10.5772/intechopen.106136*

where *ε* ¼ sup *x*∈ *D* <sup>∣</sup>*p x*ð Þ� *q x*ð Þ∣. By *<sup>I</sup>*ð Þ <sup>j</sup>*f*j*=<sup>p</sup>* denote the integral <sup>Ð</sup> <sup>∣</sup>*f x*ð Þ<sup>∣</sup> *p x*ð Þ *dx*. The both

values *ε* and *I*ð Þ j*f*j*=p* 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∪*D*2, *D*1∩*D*<sup>2</sup> ¼ ∅. For *x*∈ *D*<sup>1</sup> we get *p x*ð Þ� *q x*ð Þ� *ε*>0, and *f x*ð Þ� 0 in *D*2. If the condition *I*ð Þ¼þ j*f*j*=p* ∞ holds, then the error of the real estimation of the integral will be infinity for any *ε*> 0.

Hence, a quality of pseudorandom variables (i.e., smallness of *ε*) does not yet guaranties the smallness of the error in general, as the integral *I*ð Þ j*f*j*=p* have to be both finite and not very great in magnitude.

Suppose we simultaneously make the estimations for both *I*ð Þ j*f*j*=p* and the origin integral Eq. (2) using the same density *p x*ð Þ. Then we need to ask boundedness of the integral Ð *D* ∣*f x*ð Þ∣ *<sup>p</sup>*2ð Þ *<sup>x</sup> dx* to get the guaranteed stability of the estimation for the integral *I*ð Þ j*f*j*=p* , and so on. The qualitative comparison of simulation with both densities *p*1ð Þ *x* and *p*2ð Þ *x* 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* j*f*j*=p*<sup>1</sup> � � and *<sup>I</sup>* <sup>j</sup>*f*j*=p*<sup>2</sup> � �. On the other hand we have the Schwarz inequality:

$$\int\_{D} \frac{f(\mathbf{x})}{p(\mathbf{x})} (p(\mathbf{x}) - q(\mathbf{x})) d\mathbf{x} \le \sqrt{\int\_{D} [p(\mathbf{x}) - q(\mathbf{x})]^2 d\mathbf{x}} \sqrt{\int\_{D} \frac{f^2(\mathbf{x})}{p^2(\mathbf{x})} d\mathbf{x}}.\tag{7}$$

The sufficient condition for the estimation to be stability is that Ð *D f* 2 ð Þ *x <sup>p</sup>*2ð Þ *<sup>x</sup>* to be finite and not very great. The Schwarz inequality is reduced to the equality if and only if

$$
\lambda \frac{f(\mathbf{x})}{p(\mathbf{x})} = p(\mathbf{x}) - q(\mathbf{x}), \tag{8}
$$

where *λ* is a real number. From the above we get

$$\begin{aligned} q(\mathbf{x}) &= p(\mathbf{x}) - \lambda \frac{f(\mathbf{x})}{p(\mathbf{x})}, \\ \int\_{D} q(\mathbf{x}) d\mathbf{x} &= \int\_{D} p(\mathbf{x}) d\mathbf{x} - \lambda \int\_{D} \frac{f(\mathbf{x})}{p(\mathbf{x})} d\mathbf{x}. \end{aligned} \tag{9}$$

Therefore, the equality in Eq. (7) is reached under the necessary condition Ð *D f x*ð Þ *p x*ð Þ *dx* <sup>¼</sup> 0, when <sup>Ð</sup> *<sup>D</sup>p x*ð Þ*dx* <sup>¼</sup> <sup>Ð</sup> *<sup>D</sup>q x*ð Þ*dx* ¼ 1.

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

$$\begin{aligned} D &= [-\mathbf{1}, \mathbf{1}], \quad D\_1 = [-\mathbf{1}, \mathbf{0}], \quad D\_1 = [\mathbf{0}, \mathbf{1}], \\ f(\mathbf{x}) &= -\mathbf{1} \text{in} D\_1, \quad f(\mathbf{x}) = \mathbf{1} \text{in} D\_2, \quad p(\mathbf{x}) = p(-\mathbf{x}). \end{aligned}$$

Let us *<sup>η</sup><sup>k</sup>* be ½ � *<sup>f</sup>*ð Þ*<sup>ξ</sup> <sup>=</sup>p*ð Þ*<sup>ξ</sup> <sup>k</sup>* . Now we consider an estimation:

$$\mathcal{E}\eta\_k = \int\_D p(\mathbf{x}) \left[ \frac{f(\mathbf{x})}{p(\mathbf{x})} \right]^k d\mathbf{x} = \int\_D \frac{f^k(\mathbf{x})}{p^{k-1}(\mathbf{x})} d\mathbf{x}, \ k \ge 1. \tag{10}$$

This expectation is actually estimated by the integral: Ð *<sup>D</sup>q x*ð Þ½ � *f x*ð Þ*=p x*ð Þ *<sup>k</sup> dx*:

$$\begin{split} \int\_{D} \frac{f^{k}(\mathbf{x})}{p^{k-1}(\mathbf{x})} d\mathbf{x} - \int\_{D} q(\mathbf{x}) \frac{f^{k}(\mathbf{x})}{p^{k}(\mathbf{x})} d\mathbf{x} &= \int\_{D} [p(\mathbf{x}) - q(\mathbf{x})] \frac{f^{k}(\mathbf{x})}{p^{k}(\mathbf{x})} d\mathbf{x} \leq \\ &\leq \varepsilon \int\_{D} \frac{|f^{k}(\mathbf{x})|}{p^{k}(\mathbf{x})} d\mathbf{x}. \end{split} \tag{11}$$

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 ∣*f x*ð Þ*=p x*ð Þ∣ ≤ *M* < þ ∞ for not very great *M*. From the Schwarz inequality we have:

$$\begin{split} \int\_{D} [p(\mathbf{x}) - q(\mathbf{x})] \frac{f^{k}(\mathbf{x})}{p^{k}(\mathbf{x})} d\mathbf{x} &= \int\_{D} \frac{p(\mathbf{x}) - q(\mathbf{x})}{p^{\beta}(\mathbf{x})} \cdot \frac{f^{k}(\mathbf{x})}{p^{k-\beta}(\mathbf{x})} d\mathbf{x} \leq \\ &\leq \sqrt{\int\_{D} \frac{\left[p(\mathbf{x}) - q(\mathbf{x})\right]^{2}}{p^{2\beta}(\mathbf{x})} d\mathbf{x}} \sqrt{\int\_{D} \frac{f^{2k}(\mathbf{x})}{p^{2(k-\beta)}(\mathbf{x})} d\mathbf{x}}, \end{split} \tag{12}$$

where *β* is a real number. We get a family of proximity measures for the distribution densities:

$$\int\_{D} \frac{\left[p(\boldsymbol{\kappa}) - q(\boldsymbol{\kappa})\right]^2}{p^{2\beta}(\boldsymbol{\kappa})} d\boldsymbol{\kappa}.\tag{13}$$

For *β* ¼ 0, 5 we obtain expression

$$\chi^2(p,q) = \int\_D \frac{[p(\mathbf{x}) - q(\mathbf{x})]^2}{p^{(\mathbf{x})}d\mathbf{x}.\tag{14}$$

that is well known in the mathematical statistics. For *<sup>β</sup>* ¼ �0, 5 we get <sup>Ð</sup> *<sup>D</sup>*½ � *p x*ð Þ� *q x*ð Þ <sup>2</sup> *p x*ð Þ*dx* and have the obvious inequalities

$$\begin{aligned} \int\_{D} \left[ p(\mathbf{x}) - q(\mathbf{x}) \right]^2 p(\mathbf{x}) d\mathbf{x} &\leq \sup\_{\mathbf{x} \in D} p(\mathbf{x}) \cdot \int\_{D} \left[ p(\mathbf{x}) - q(\mathbf{x}) \right]^2 d\mathbf{x}, \\\int\_{D} \left[ p(\mathbf{x}) - q(\mathbf{x}) \right]^2 p(\mathbf{x}) d\mathbf{x} &\leq \sup\_{\mathbf{x} \in D} \left[ p(\mathbf{x}) - q(\mathbf{x}) \right]^2 \int\_{D} p(\mathbf{x}) d\mathbf{x} = \sup\_{\mathbf{x} \in D} \left[ p(\mathbf{x}) - q(\mathbf{x}) \right]^2. \end{aligned}$$

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

$$\frac{p-q}{p^{\beta}} = \lambda \frac{f^k}{p^{k-\beta}}, \quad p-q = \lambda \frac{f^k p^{\beta}}{p^{k-\beta}},\tag{15}$$

i.e., the necessary condition is Ð *D f k* ð Þ *x pk*�2*<sup>β</sup>* ð Þ *<sup>x</sup> dx* <sup>¼</sup> <sup>0</sup>*:* This is realized in the Example 2. Let us remark that for *k* ¼ 1 and *β* ¼ 1 we have

$$\begin{split} \int\_{D} \frac{p(\mathbf{x}) - q(\mathbf{x})}{p(\mathbf{x})} f(\mathbf{x}) d\mathbf{x} &= \int\_{D} \frac{p - q}{\sqrt{p}} \cdot \frac{f}{\sqrt{p}} \leq \\ &\leq \sqrt{\int\_{D} \frac{(p - q)^{2}}{p} d\mathbf{x}} \sqrt{\int\_{D} \frac{f^{2}}{p} d\mathbf{x}} = \sqrt{\chi^{2}(p, q)} \cdot \mathcal{E} \eta\_{2}. \end{split} \tag{16}$$

We assume that integrals Ð *<sup>D</sup>*½ � *f x*ð Þþ *εi*ð Þ *x dx* are known and the subintegral functions *f x*ð Þþ *εi*ð Þ *x* >0 close to a function *f x*ð Þ≥0. Suppose also

$$\begin{aligned} 1 - \delta\_1(\varepsilon\_i) &< \frac{f(\mathbf{x})}{f(\mathbf{x}) + e\_i(\mathbf{x})} < 1 + \delta\_2(\varepsilon\_i),\\ 1(f) - \delta\_3(\varepsilon\_i) &\le & (f + e\_i) \le f(f) + \delta\_4(e\_i), \quad \delta\_j \to 0, \quad \text{as} \quad e\_i \to 0, \ j = 1, 2, 3, 4. \end{aligned} \tag{17}$$

$$\int\_D q(\mathbf{x}) \frac{f(\mathbf{x})}{p(\mathbf{x})} d\mathbf{x} = \int\_D q(\mathbf{x}) \frac{f(\mathbf{x}) f(f + e\_i)}{f(\mathbf{x}) + e\_i(\mathbf{x})} d\mathbf{x} = \mathcal{E}\hat{\eta},$$

where the random value ^*η* has a form:

$$
\hat{\eta} = \frac{f\left(\hat{\xi}\right)f(f+\varepsilon\_i)}{f\left(\hat{\xi}\right) + \varepsilon\_i\left(\hat{\xi}\right)}.\tag{18}
$$

Here ^*<sup>ξ</sup>* is distributed with the density *q x*ð Þ. Keeping the above factors in mind we get the following:

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


Let we calculate Ð *D* ∣*f x*ð Þ∣ *p x*ð Þ *dx* using the density

$$p\_1(\mathbf{x}) = \frac{|f(\mathbf{x})|}{I(|f|/p)p(\mathbf{x})},$$

then the estimation variance equals to zero.

Suppose we calculate the integral Ð *<sup>D</sup>f x*ð Þ*dx* with the density *p*1ð Þ *x* . In this case it would be interesting to know both the values Ð *D* ∣*f x*ð Þ∣ *p*1 *dx* and Ð *D f* 2 *p*2 1 *dx*.

**Proposition 1.** *I* j*f*j*=p*<sup>1</sup> � � <sup>¼</sup> *<sup>I</sup>*ð Þ <sup>j</sup>*f*j*=<sup>p</sup> .* **Proposition 2.** *I f* <sup>2</sup> *=p*<sup>2</sup> 1 � � <sup>¼</sup> *<sup>I</sup>* 2 ð Þ <sup>j</sup>*f*j*=<sup>p</sup> I p*<sup>2</sup> ð Þ*.* Now we consider the density

$$p\_2(\mathbf{x}) = \frac{f^2(\mathbf{x})}{p(\mathbf{x})I(f^2/p)}.$$

Using the density above for the estimation of the integral Ð *D f* 2 *p dx* we obtain the estimation variance equals to zero.

Further we estimate the integral Ð *<sup>D</sup>f x*ð Þ*dx* with the density *p*2: ^*η* ¼ *f ξ*<sup>2</sup> ð Þ*=p ξ*<sup>2</sup> ð Þ, *ξ*<sup>2</sup> is distributed with *p*2ð Þ *x* . Suppose *η* ¼ *f*ð Þ*ξ =p*ð Þ*ξ* , *ξ* is distributed with *p x*ð Þ; then *I f* <sup>2</sup> *=p*<sup>2</sup> � � <sup>¼</sup> *I f* <sup>2</sup> *=p* � �.

**Proposition 3.** *varη* ¼ *var*^*η.*

$$\textbf{Proposition 4.}\text{ }[\underset{D}{\underset{D}{\underset{p^{(x)}}{p^{(x)}}}}]{\textbf{1}.}\text{ }[\underset{D}{\underset{p^{(x)}}{\right]}}\text{ }\leq(vol\,D)\cdot\underset{D}{\underset{p^{(x)}}{\rightleftharpoons}}\text{ }\forall\text{x},\text{where }vol\,D\text{ is the volume of the}$$

domain *D*.

**Proposition 5.** If *volD* <sup>¼</sup> 1, *<sup>p</sup>*3ð Þ¼ *<sup>x</sup> <sup>f</sup>* 2 ð Þ *x I f* <sup>2</sup> ð Þ, *<sup>η</sup>*<sup>3</sup> <sup>¼</sup> *<sup>f</sup> <sup>ξ</sup>*<sup>3</sup> ð Þ *<sup>p</sup> <sup>ξ</sup>*<sup>3</sup> ð Þ , where *<sup>ξ</sup>*<sup>3</sup> is a random variable distributed with the density of *<sup>p</sup>*3ð Þ *<sup>x</sup>* ; *<sup>η</sup>*<sup>4</sup> <sup>¼</sup> *<sup>f</sup> <sup>ξ</sup>*<sup>4</sup> ð Þ *<sup>p</sup> <sup>ξ</sup>*<sup>4</sup> ð Þ , where *<sup>ξ</sup>*<sup>4</sup> is a random variable distributed with the density of *p x*ð Þ� 1, then *varη*<sup>3</sup> ¼ *varη*4.

In actual practice normalization constants are usually unknown for both *p*1ð Þ *x* and *p*2ð Þ *x* . But using densities close to them we can get the approximate equalities in the Prepositions 1, 2, 3.

Now we consider

$$I(f) = \int\_{[0,1]^{10}} \mathbf{x}\_1 \mathbf{x}\_2 \dots \mathbf{x}\_{10} \, d\mathbf{x}\_1 \dots d\mathbf{x}\_{10},\tag{19}$$

where the integration domain *<sup>D</sup>* <sup>¼</sup> ½ � 0, 1 <sup>10</sup> is the 10-dimensional unit cube, the subintegral function *f x*ð Þ is equals to *x*1*x*<sup>2</sup> … *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 *η* ¼ *f*ð Þ*ξ =p*ð Þ*ξ* , where *ξ* is distributed with the density *p*ð Þ*ξ* . In fact we realize the discrete set of numbers *ξi*, *i* ¼ 1, … , *N*, which we can consider to be realizations of some distribution *qN*ð Þ *x* .

In all our numerical computations we use the pseudorandom number generator: *generator type* mt19937 [3]. For *s* ¼ 10, *p x*ð Þ� 1 we have

$$\begin{aligned} I(f) &= (\mathbf{1}/2)^{10} \approx \mathbf{9}, \mathbf{7}656 \cdot \mathbf{10}^{-4}, \quad I(f/p) = I(f), \\ I(f^2/p) &= I(f^2/p^2) = (\mathbf{1}/3)^{10} \approx \mathbf{1}, 6935 \cdot \mathbf{10}^{-5}, \\ \sigma &= \sqrt{(\mathbf{1}/3)^{10} - (\mathbf{1}/2)^{20}} \approx \mathbf{3}, 998 \cdot \mathbf{10}^{-3}. \end{aligned}$$

**Table 1** shows the empirical estimations ^*I* and *σ*^ for *I* and *σ*, respectively. Taking *p x*ð Þ¼*<sup>i</sup>* <sup>3</sup>*x*<sup>2</sup> *<sup>i</sup>* over each coordinate we get

$$I(\mathfrak{f}^2/p) = \left(\mathbf{1}/\mathfrak{3}\right)^{10}, \quad I(\mathfrak{f}/p) = \infty, \quad I(\mathfrak{f}^2/p^2) = \infty,$$

$$\sigma = \sqrt{\left(\mathbf{1}/\mathfrak{3}\right)^{10} - \left(\mathbf{1}/\mathfrak{2}\right)^{20}}.$$

The value *σ* is the same as one for *p x*ð Þ� 1.


**Table 1.**

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

*Stability of Algorithms in Statistical Modeling DOI: http://dx.doi.org/10.5772/intechopen.106136*


**Table 2.**

*The results of numerical calculations for the integral Eq. (19) with the density p x*ð Þ¼ *<sup>i</sup>* <sup>3</sup>*x*<sup>2</sup> *i* .

For *N* ¼ 100000000 the computer code outputs an error because of machine zero divide. The reason of this event is *<sup>η</sup>* <sup>¼</sup> <sup>Q</sup><sup>10</sup> <sup>1</sup> 1*=* 3 ffiffiffi *α* p<sup>3</sup> *i* � � �, where *α<sup>i</sup>* are the pseudorandom numbers with the uniform density in 0, 1 ð Þ [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 f*ð Þ is successfully estimated, but the empirical estimations of *σ*^ are sufficiently different from the theoretical value *σ*. This result is explained by the following: <sup>E</sup>*η*<sup>4</sup> <sup>¼</sup> <sup>∞</sup>, *I f*ð Þ¼ *<sup>=</sup><sup>p</sup>* <sup>∞</sup>, *I f* <sup>2</sup> *<sup>=</sup>p*<sup>2</sup> � � <sup>¼</sup> <sup>∞</sup>*:* Taking *p x*ð Þ¼*<sup>i</sup>* 2, 6*x*1,6 *<sup>i</sup>* we obtain *<sup>σ</sup>* <sup>≈</sup>1, 223 � <sup>10</sup>�3, *I f*ð Þ *<sup>=</sup><sup>p</sup>* <sup>≈</sup> 0, 67556, the finite value of <sup>E</sup>*η*4, and *I f* <sup>2</sup> *<sup>=</sup>p*<sup>2</sup> � � <sup>¼</sup> <sup>∞</sup>. Although the last estimation is infinite, but **Table 3** shows that both values *I f*ð Þ and *σ*^ are successfully calculated. The value of *σ*^ is very close to *σ*.
