How to calculate joint distribution by integrating over all possible values of model parameters and observations

Question

I was wondering how to calculate the following joint distribution by assuming that $x_i$'s are continuous observations from normal distirbution $N(\mu, \sigma)$ with mean $\mu$ and variance $\sigma$ and $p(\mu, \sigma)$ is the joint distribution of $\mu$ and $\sigma$: $$ p(x_1,x_2,..., x_n, x_{n+1}) = \iint {\displaystyle \prod_{j=1}^{n+1} N(\mu,\sigma)|_{x_j}} . p(\mu, \sigma) \,d\mu \, d\sigma $$ I know that this can be accomplished by solving the intergration numerically (e.g. with rectangular integration), however, the thing that puzzles me is that the products of $N(\mu,\sigma)|_{x_j}$ would be always zero since the probability of the point $x_i$ in a continuous distribution $N(\mu,\sigma)$ is zero (or I'm wrong?).

In addition, what is the best method to calculate (or estimate) this distribution in python?

Answer 1

Ok, this looks like deep learning related stuff, maybe an exercise.

$N(\mu,\sigma)|_{x_j}$ is not the probability of a point $x_j$, this is the density of point $x_j$ given the parameters $\mu, \sigma$. The total integral is not the probability either, it is the likelihood of sorts or a density function. The whole thing is in fact never a zero, because the density $N$ of normal distribution is never zero, it can be very small of course.