Estimating covariance of the difference of directional distributions derived from Gaussian mixtures Given Gaussian mixtures $X_1, X_2 \in \mathbb{R}^p$ defined as $$P(X_i = x) = \sum_s \omega^{(s)}_i \mathcal{N}(x; \mu^{(s)}_i, \Sigma_i)$$ where the superscript $(s)$ indexes the $s$-th component of the mixture distribution $X_i$, the covariance of all mixture components are identical, and $\sum_s \omega^{(s)}_i = 1$ to keep the resulting density normalized. Define new random variables $Z_i$ having support only on the unit sphere in $\mathbb{R}^p$ by marginalizing $X_i$ over the distance from the origin,
\begin{equation}
P(Z_i = z) =
\begin{cases}
\int^\infty_0 P(X_i = t z)\,dt & \|z\| = 1 \\
0 & \|z\| \neq 1 .
\end{cases}
\end{equation}
I am looking to find computationally efficient methods for estimating $\Gamma = \mathbf{E}[(Z_1 - Z_2)(Z_1 - Z_2)^T]$.
I have thought of two approaches so far:


*

*Crudest possible approach: estimate a Gaussian $Y_i$ that best approximates $Z_i$ (by performing importance sampling with $X_i$ as the proposal distribution to compute $Z_i$'s mean and covariance), then use $\mathbf{E}[(Y_1 - Y_2)(Y_1 - Y_2)^T]$ as the estimate of $\Gamma$.

*Perform importance sampling directly on $Z_1 - Z_2$ by drawing samples $x_1, x_2$ from $X_1, X_2$, computing the corresponding values of $z_1, z_2$, and the importance weight as $P(Z_1 = z_1, Z_2 = z_2)$.
 A: Approach (1) might be feasible if the magnitude of the covariances are small compared to the (square of) the magnitude of the mean vectors.  If this is the case, then the projection of the $X_i$ onto the unit sphere will subtend a small angle (this ignores some of the data in the tails), and you can ignore the curvature of the sphere. 
For approach (2) I'm not sure that importance sampling is really required;
it is easy to draw a random number from a Gaussian mixture model:


*

*Pick a kernel using the prior over kernels, i.e. from the $\omega_i^{s}$,

*Draw a Gaussian number from that kernel, which is easy enough if you
find a library to find the square-root of the covariance matrix.


Then each sample has a weight of $1$.
If your problem is that you don't want to find the square-root of each of the covariance matrices, I'd use a sampling distribution that is centered on the given kernel, spherically symmetric in $\Re^p$, and whose scale is set by requiring that the determinant of the symmetric sampling distribution is equal to that of the true distribution.  Then the importance weight of each sample is just the ratio of the true distribution's pdf to that of this symmetric sampling distribution. 
In the end it might be useful to think of the computation in terms of a table:

       k_1    k_2    k_3    ...   | 
==================================+
       M_11   M_12   M_13   ...   |  K_1
       M_21   ...                 |  K_2
       ...                        |  ...


Where $M_{st}$ represents the value of the moment $E[ (Z_1-Z_2)^2]$ for the s'th  $X_1$ kernel (i.e. pick the index $s$ in distribution 1) and $X_2$ is drawn from
the t'th kernel of the $X_2$ distribution.  Then the overall moment is
$$
\sum_{st} \omega_1^{(s)} \omega_2^{(t)} M_{st}
$$
A: A potentially useful result for circularly symmetric Gaussian kernels in 3-D:
First define:
$$ 
erf(x)=\int_{0}^{x} \frac{2}{\sqrt{\pi}} e^{-t^2} dt 
$$
$$
Q(x)=\int_{0}^{x} \frac{2}{\sqrt{\pi}} t^2 e^{t^2} dt = \frac{1}{2} erf(x)-\frac{x}{\sqrt{\pi}} e^{-x^2} 
$$
The following integral can be evaluated in terms of these
$$
I(\bar{x}, \sigma) =\int_{0}^{\infty} x^2 exp \left[ -\frac{1}{2} \frac{ (x-\bar{x})^2}{\sigma^2} \right] dx
$$
via
$$ 
u=\frac{1}{\sqrt{2}} \frac{ (x-\bar{x})}{\sigma}\
$$
$$
I(\bar{x}, \sigma)=\sqrt{\frac{\pi}{2}} \sigma  \left[ 2 \sigma^2 \left( \frac{1}{2} -Q( \frac{-\bar{x}}{\sqrt{2} \sigma} ) \right) + 2 \sqrt{ \frac{2}{\pi}} \sigma \bar{x} exp \left[ \frac{-\bar{x}}{\sqrt{2} \sigma} \right] +\bar{x}^2 \left( 1-erf(\frac{-\bar{x}}{\sqrt{2}\sigma} ) \right) \right]
$$
If we have 3-D isotropic Gaussian distribution:
$$
p(\vec{u} \vert \vec{\mu} , \sigma^2 ) = (2 \pi \sigma^2 )^{-3/2} exp \left[ -\frac{1}{2} \frac{ (\vec{u}-\vec{\mu})^2} {\sigma^2} \right]
$$
And define the probability density as a function of direction, here
expressed in terms of the unit vector $\hat{v}$
$$
p(\hat{v} \vert \vec{\mu} , \sigma) = \int_0^{\infty} r^2 dr (2 \pi \sigma^2)^{-3/2} exp \left[ -\frac{1}{2} \frac{ (r \hat{v} -\vec{\mu} )^2}{\sigma^2} \right] 
$$
And by completing the square in the exponent you get
$$
p(\hat{v} \vert \vec{\mu} , \sigma ) = (2 \pi \sigma^2)^{-3/2} exp \left[ 
- \frac{1}{2} \frac{ \vec{\mu} \cdot \vec{\mu} }{\sigma^2}
+ \frac{ (\hat{v} \cdot \vec{\mu} )^2 }{ \sigma^2 }
\right] I( \hat{v} \cdot \vec{\mu} , \sigma )
$$
As mentioned above, this is derived for the 3-D case, but hopefully it gives an idea how to go about it for the 2-D case.
Thus, if you know (or estimate) the mean and the variance for each kernel, you can derive the distribution for the direction vectors, and thus, compute 
it's moments.
