Expectation of function of random variables (example Euclidean distance)

Question

How the expectation value of a nonlinear transform of normal distributed variables can be approximated by Mathematica in case there is no analytical solution?

Let's take the Euclidian distance of 2 points in $\mathbb{R^3}$ as an example. The variables $p0,...,p5$ shall be the expectation values of the point coordinates. The distributions are independent and identical with equal variance $s^2$.

Below is the code where an approximated solution in dependence of $p0,...,p5$ and $s$ shall be found:

Expectation[
 Sqrt[(x0 - x1)^2 + (x2 - x3)^2 + (x4 - x5)^2],
 {
  x0 \[Distributed] NormalDistribution[p0, s^2],
  x1 \[Distributed] NormalDistribution[p1, s^2],
  x2 \[Distributed] NormalDistribution[p2, s^2],
  x3 \[Distributed] NormalDistribution[p3, s^2],
  x4 \[Distributed] NormalDistribution[p4, s^2],
  x5 \[Distributed] NormalDistribution[p5, s^2]
  }]

Answer 1

From @user64494 's comments and the mentioned [Wiki page][1] the expectation of

$$\sqrt{(x_0-x_1)^2+(x_2-x_3)^2+(x_4-x_5)^2}$$

will be

$$\sqrt{\pi } \sigma L_{\frac{1}{2}}^{\frac{1}{2}}\left(-\frac{L^2}{4 \sigma ^2}\right)$$

or in Wolfram Language

σ Sqrt[π]*LaguerreL[1/2, 1/2, -L^2/(4 σ^2)]

where I've substituted the more commonly used $\sigma$ for $s$ and

L = Sqrt[(p0 - p1)^2 + (p2 - p3)^2 + (p4 - p5)^2]

While this is an exact result the OP has the concern that few programming languages have adequate approximations (if they have an approximation at all) of the LaguerreL function.

Fortunately, there is a simplification in this case.

σ Sqrt[π]*LaguerreL[1/2, 1/2, -L^2/(4 σ^2)] // FunctionExpand // FullSimplify

results in

$$\left(\frac{2 \sigma ^2}{L}+L\right) \text{erf}\left(\frac{L}{2 \sigma }\right)+\frac{2 \sigma e^{-\frac{L^2}{4 \sigma ^2}}}{\sqrt{\pi }}$$

or

(2 E^(-(L^2/(4 σ^2))) σ)/Sqrt[π] + (L + (2 σ^2)/L) Erf[L/(2 σ)]

Most languages (even Excel) have an Erf function (or an equivalent) readily available.

(Again, this is not an approximation as requested. If the OP can present an example where there isn't an exact solution, then approximation approach can be better targeted. My reasoning for this is that such approximations depend on the functions being approximated and the range of values for which the approximation is desired to work. Being specific for a needed approximation would likely get better answers.)

Answer 2

There are simpler solutions for special cases.

Definitions

$x_i=\mathcal{N}(p_i,\sigma^2)_{i=1...6}$ are 6 independent normal distributed random variates with expectations $p_i$ and variance $\sigma^2$ (this replaces $s^2$ from the original post).

The random variables are transformed by the Euclidian distance function $$L(x_0,x_1,x_2,x_3,x_4,x_5)=\sqrt{(x_0 - x_1)^2 + (x_2 - x_3)^2 + (x_4 - x_5)^2}$$ The distance of the expected values is $$L(p_0,p_1,p_2,p_3,p_4,p_5)=L_p=\sqrt{(p_0 - p_1)^2 + (p_2 - p_3)^2 + (p_4 - p_5)^2}$$

The expected value of $L$ is $\mathbb{E}[L]$.

Results

For following 2 cases simpler expressions of approximated expectations $\mathbb{E_1}[L], \mathbb{E_2}[L] $ are:

$$\frac{\sigma}{L_p} \approx 0 \hspace{2mm} \rightarrow \hspace{2mm} \mathbb{E_1}[L] =\sqrt{L_p^2+6\sigma^2} \approx \mathbb{E}[L] $$

$$\frac{L_p}{\sigma} \approx 0 \hspace{2mm} \rightarrow \hspace{2mm} \mathbb{E_2}[L]=\frac{4\sigma}{\sqrt{\pi}}\approx \mathbb{E}[L] $$

Derivations

The distribution $$L'(x_i)=\left(\frac{x_0 - x_1}{\sqrt{2}\sigma}\right)^2+\left(\frac{x_2 - x_3}{\sqrt{2}\sigma}\right)^2+\left(\frac{x_4 - x_5}{\sqrt{2}\sigma}\right)^2$$ is a noncentral chi square distribution with $k=3$ degrees of freedom and the noncentrality parameter is

$$\lambda=L'(p_i)=\left(\frac{p_0 - p_1}{\sqrt{2}\sigma}\right)^2+\left(\frac{p_2 - p_3}{\sqrt{2}\sigma}\right)^2+\left(\frac{p_4 - p_5}{\sqrt{2}\sigma}\right)^2=\left(\frac{L_p}{\sqrt{2}\sigma}\right)^2$$

The expectation value of a noncentral chi square distribution is

$$\mathbb{E}[L']=\lambda+k$$

If $\frac{\sigma}{L_p}\approx 0 $ then $L' \approx \lambda$ and using the approximation $$\frac{1}{n}\sum_{i=1}^n \sqrt{\lambda+\delta_i}\approx \sqrt{\frac{1}{n}\sum_{i=1}^n (\lambda+\delta_i)} \hspace{3mm}\text{for}\hspace{3mm} |\delta_i|\ll \lambda $$

$\mathbb{E}[L]$ is approximated by

$$\frac{\mathbb{E_1}[L]}{\sqrt{2}\sigma}=\mathbb{E_1}\left[\sqrt{L'}\right]\approx \sqrt{\mathbb{E}[L']}$$ and it follows $$\mathbb{E_1}[L] =\sqrt{L_p^2+6\sigma^2}$$

If $\frac{L_p}{\sigma}\approx 0 $ then $\lambda\approx 0$ and $\mathbb{E}[L]$ is approximated by the expectation of a chi distribution

$$\frac{\mathbb{E_2}[L]}{\sqrt{2}\sigma}\approx\sqrt{2}\frac{\Gamma\left(\frac{k+1}{2}\right)}{\Gamma\left(\frac{k}{2}\right)}=2\sqrt{\frac{2}{\pi}}$$ It follows $$\mathbb{E_2}[L]=\frac{4\sigma}{\sqrt{\pi}}$$

It holds $\mathbb{E_1}[L] \ge \mathbb{E}[L] \ge \mathbb{E_2}[L] $.

Examples

If the values of $L_p$ and $\sigma$ are in the proper ratio then the approximations are quite close to the exact solutions. This is seen below where the approximate and exact solutions are compared. The example coordinates from another post $p_i=(1,0,2,0,3,0)$ are used, $L_p=\sqrt{14}\approx 3.74$.

• $\sigma=0.02, \frac{\sigma}{L_p}\approx 0$

$\small \mathbb{E_1}[L]=3.7419781$

$\small\mathbb{E}[L]=3.7418712$

• $\sigma=200, \frac{L_p}{\sigma}\approx 0$

$\small\mathbb{E_2}[L]=451.35167$

$\small\mathbb{E}[L]=451.36483$

• $\sigma=2 \approx L_p $

$\small\mathbb{E_1}[L]=6.1644140$

$\small\mathbb{E_2}[L]=4.5135167$

$\small\mathbb{E}[L]=5.7275959$

Motivation

Why do we need these approximations? The exact solutions require calculation of the mean of a noncentral chi distribution that practically means calculation of a generalized Laguerre function that is related to Kummer's confluent hypergeometric function. In Mathematica the functions $Hypergeometric1F1[a,b,z]$ with $a=-\frac{1}{2}, b=\frac{3}{2}$ or $LaguerreL[n,k,z]$ with $n=k=\frac{1}{2}$ can be used, $z$ has negative real part. There might be problems if you try to translate your code into another languange. Such special functions are not available in all software packages or supported only for $a,b,n,k\in \mathbb{N}$. Special treatment regarding numerical stability is needed if you write your own routine. So the given approximations may help in some cases.