2
$\begingroup$

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]
  }]
$\endgroup$
8
  • $\begingroup$ Is p0, p1,... and s not known? $\endgroup$ – MikeY Jan 16 '20 at 18:31
  • $\begingroup$ These are variables. $\endgroup$ – granular bastard Jan 16 '20 at 18:44
  • $\begingroup$ This is rather math than Mathematica. Because the difference of two independent normal distributions is a normal distribution too, you deal with the mean of a noncentalchisquare distribution( see en.wikipedia.org/wiki/Noncentral_chi_distribution and reference.wolfram.com/language/ref/… ). $\endgroup$ – user64494 Jan 16 '20 at 19:20
  • 1
    $\begingroup$ Your claim "The noncentral cchi distribution has no closed solution for the expectation value" is not exact Its mean is expressed it terms of the Laguerre function (see the above linked Wiki article) so this is an analytic expression (see en.wikipedia.org/wiki/Closed-form_expression). MMA deals with Laguerre functions. Don't hesitate to ask for further explanation in need $\endgroup$ – user64494 Jan 16 '20 at 19:47
  • 1
    $\begingroup$ No, L[1/2,1/2,x] is not expressed through other special and elementary functions. However, we can built its plot, calculate LaguerreL[1/2, 1/2, 2.3] in MMA and so on. $\endgroup$ – user64494 Jan 16 '20 at 20:38
9
$\begingroup$

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.)

$\endgroup$
6
  • $\begingroup$ Please look here stats.stackexchange.com/questions/445185/… $\endgroup$ – granular bastard Jan 16 '20 at 23:59
  • $\begingroup$ Can Mathematica help if the solution is not known? $\endgroup$ – granular bastard Jan 20 '20 at 9:04
  • $\begingroup$ Certainly. Finding the mean after using Series on the random variable for an approximation (or if one is lucky that can find the exact mean). Or if there are few enough parameters, use numerical integration for a wide range of parameters and then see if there is a function of the parameters that will predict the mean. And setting a specific value for each parameter and then taking the mean of lots of random samples. The latter is easily programmed in most languages. $\endgroup$ – JimB Jan 20 '20 at 17:57
  • $\begingroup$ If $\sigma=0.001$ and L is in the range [0,1] then z is in the range [0,-250000]. There is no practical series expansion for this range. One needs an expansion in dependence of $u$, however this is not calculated: $Series[LaguerreL[1/2, 1/2, z], {z, u, 1}]=LaguerreL[1/2,1/2,u]-LaguerreL[-(1/2),3/2,u] (z-u)+O[z-u]^2$ $\endgroup$ – granular bastard Jan 20 '20 at 21:51
  • $\begingroup$ Not true. I just stumbled across this: FunctionExpand[LaguerreL[1/2, 1/2, -z]] // FullSimplify results in $\frac{(2 z+1) \text{erf}\left(\sqrt{z}\right)}{\sqrt{\pi } \sqrt{z}}+\frac{2 e^{-z}}{\pi }$. I'll add that into my answer. $\endgroup$ – JimB Jan 21 '20 at 16:27
-3
$\begingroup$

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.

$\endgroup$
9
  • $\begingroup$ The statement $\frac{\sigma}{L} \approx 0 $ makes no sense: its LHS is a random variable, whereas its RHS is a number. Also the relation $ \approx $ is not defined. $\endgroup$ – user64494 Jan 19 '20 at 9:29
  • $\begingroup$ The statement after edit $\frac{\sigma}{L(p_i)} \approx 0 $ makes no sense because of the same reasons. $\endgroup$ – user64494 Jan 19 '20 at 10:44
  • $\begingroup$ This answer does not deal with Mathematica at all. Next, it looks like a fake: several edits in the same formula, no definitions, no calculations, and no references. $\endgroup$ – user64494 Jan 19 '20 at 11:09
  • $\begingroup$ The answer was refurbished and includes now proofs and motivation. $\endgroup$ – granular bastard Jan 19 '20 at 15:07
  • $\begingroup$ Up to en.wikipedia.org/wiki/Noncentral_chi_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 not a noncentral chi distribution: there is no square root. Also there is $x_i$ in the LHS, but there is no $x_i$ in the RHS. $\endgroup$ – user64494 Jan 19 '20 at 15:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.