# Expectation of the Gaussian likelihood

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

• The term "likelihood" is not used correctly in this question. It should not appear since variability over the parameter values is not of a concern. Commented Dec 10, 2023 at 8:53

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

• Hi @Cryo, thanks for your reply and well worked solution, really appreciate the feedback. Re-writing it in quadratic form is very clever. However, can I please clarify your solution? It doesn't seem to match up with my simulation. For simplicity, take u2=0 and u1=0 so that the aB and gamma terms drop out. Eg. in python, 1000 draws from a Z ~ N(0,2) and then calculating the expected likelihood mean(L(Z | u1=0, sd1=1)), I get approx 0.184. With the above formula I get 0.226. As fair as I can tell, something seems to be off. I did the calculation empirically in Python. Commented Dec 10, 2023 at 0:10
• @reynolds.brian yep, my bad - forgot square in the definition of $\alpha$. Fixed now
– Cryo
Commented Dec 10, 2023 at 8:42
• Amazing, thankyou! Really appreciated. I'd like to make sure you're credited properly, I'll put my email and LinkedIn on my public profile on my for a few days if you would like to get in contact. Otherwise, thanks once again. Commented Dec 10, 2023 at 10:51