Expectation of the Gaussian likelihood

Question

I'm working on a challenging machine learning problem, where I need to find the expectation of the likelihood of one Gaussian, given the parameters of another. Apologies if any of the notation is rough, this is new territory for me.

More precisely, given two independent Gaussian distributions.

$$X \sim \mathcal{N}(\mu_1, \sigma_1^2) $$

$$Z \sim \mathcal{N}(\mu_2, \sigma_2^2)$$

where the values of $\mu_1$, $\sigma_1$, $\mu_2$ and $\sigma_2$ are known. Suppose I were to randomly draw values from $Z$ and evaluate the likelihood using the pdf of $X,$ $p(X)$ with $\mu_1, \sigma_1$, what is the expected likelihood? To clarify, not the expected value of $Z$ but the expectation of the likelihood function itself.

I've managed to get the expected likelihood $E[p(X)]$, which holds for $Z$ when $\mu_1=\mu_2$ and $\sigma_1=\sigma_2$, via the following:

Definition of the pdf of $X$ is:

$$ p(X) = \int_{-\infty}^\infty \frac{1}{\sqrt{2\pi}\sigma_1} e^{-\frac{(x-\mu_1)^2}{2\sigma_1^2}} \, dx$$

and the definition of the expectation is:

$$E[x] = \int_{-\infty}^\infty x \cdot f(x) \ dx$$ $x = p(X)$ so $$E[p(X)] = \int_{-\infty}^\infty p(X) \cdot p(X) \ dx$$

$$E[p(X)] = \int_{-\infty}^\infty \left(\frac{1}{\sqrt{2\pi}\sigma_1} e^{-\frac{(X-\mu_1)^2}{2\sigma_1^2}}\right) \cdot \left(\frac{1}{\sqrt{2\pi}\sigma_1} e^{-\frac{(X-\mu_1)^2}{2\sigma_1^2}}\right) \ dx$$

rearranging a little gives:

$$E[p(X)] = \frac{1}{\sigma_1^2 \cdot 2 \pi} \cdot \int_{-\infty}^\infty \exp\left(-\frac{(x - \mu_1)^2}{\sigma_1^2}\right) \ dx $$

and finally the following result, which I've verified empirically via simulation:

$$ E[p(X)] = \frac{1}{\sigma_1^2 \cdot 2 \pi} \cdot {\sigma_1 \sqrt{\pi}} $$

But, how do I calculate the likelihood when integrating over $Z$, with $\mu_1$ and $\sigma_1$ from $X$?

I tried defining $x = \frac{(z - \mu_2)\sigma_1}{\sigma_2}$ and using the change of variables approach, substituting and integrating over $dz.$ However, that doesn't work, I'm definitely doing something wrong.

How should I proceed from here? I'm lost and any help would be much appreciated.

Thanks

Answer 1

Firstly, I don't think your definition of $p\left(x\right)$ is correct. You need to give it without the integral, but this is a minor thing.

You mostly did the right things after that, I think. So what you are trying to evaluate is:

$$ \begin{align} E_z\left[p_x\left(z\right)\right]&=\frac{1}{\sigma_2\sqrt{2\pi}}\int dz \exp\left(-\frac{\left(z-\mu_2\right)}{2\sigma_2^2}\right)\,p_x\left(z\right)=\\ &=\frac{1}{\sigma_2\sqrt{2\pi}}\int dz \exp\left(-\frac{\left(z-\mu_2\right)^2}{2\sigma_2^2}\right)\,\left\{\frac{1}{\sigma_1\sqrt{2\pi}}\,\exp\left(-\frac{\left(z-\mu_1\right)^2}{2\sigma_1^2}\right)\right\} \end{align} $$

You can re-write this as:

$$ \begin{align} \alpha^2&=\frac{1}{2\sigma_2^2}+\frac{1}{2\sigma_1^2}\\ \alpha\beta&=\frac{\mu_2}{2\sigma_2^2}+\frac{\mu_1}{2\sigma_1^2}\\ \gamma&=\frac{\mu_2^2}{2\sigma_2^2}+\frac{\mu_1^2}{2\sigma_1^2}\\ E_z\left[p_x\left(z\right)\right]&=\frac{1}{\sigma_1\sigma_2\cdot 2\pi}\int dz \exp\left(-\left(\alpha^2 z^2 - 2\alpha\beta z +\gamma\right)\right) \end{align} $$

You can then proceed to complete the square: $$ \begin{align} E_z\left[p_x\left(z\right)\right]&=\frac{1}{\sigma_1\sigma_2\cdot 2\pi}\cdot\exp\left(-\gamma+\beta^2\right)\int dz \exp\left(-\left(\alpha z - \beta\right)^2\right)=\\ &=\frac{1}{\sigma_1\sigma_2\cdot 2\pi}\cdot\exp\left(-\gamma+\beta^2\right)\cdot\frac{1}{\alpha}\int d\zeta \exp\left(-\zeta^2\right)\end{align} $$

Where the last step is a simple variable substitution. The integral is standard, the rest of the variables are known

Thus:

$$ \begin{align} E_z\left[p_x\left(z\right)\right]&=\frac{1}{2\sigma_1\sigma_2\alpha\cdot \sqrt{\pi}}\cdot\exp\left(-\gamma+\beta^2\right) \end{align} $$

---

import numpy as np
import scipy.stats as sp_st

mu1 = 0.2
mu2 = 0.3
#
sig1 = 0.4
sig2 = 0.5
#
sample_size = 10000000

###### sampled value
sample2 = sp_st.norm(loc=mu2, scale=sig2).rvs(sample_size)
samp_expectation2_p1 = np.mean( 
    sp_st.norm(loc=mu1, scale=sig1).pdf(sample2) 
)
print(f'Sampled expectation2_p1: {samp_expectation2_p1:.3e}')

##### analystic value
alpha2 = 1/(2*(sig1**2)) + 1/(2*(sig2**2))
alpha = np.sqrt(alpha2)
alpha_beta = mu1/(2*(sig1**2)) + mu2/(2*(sig2**2))
beta = alpha_beta/alpha
gamma = (mu1**2)/(2*(sig1**2)) + (mu2**2)/(2*(sig2**2))
#
analytic_expectation2_p1 = (1/(2*sig1*sig2*alpha*np.sqrt(np.pi)))*\
    np.exp(-gamma+beta**2)
print(f'Analytic expectation2_p1: {analytic_expectation2_p1:.3e}')

Output:

Sampled expectation2_p1: 6.156e-01
Analytic expectation2_p1: 6.155e-01