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

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?

$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.
• Thanks. Do you have any suggestion for the implementation in python? I've used dblquad for the integration but it takes a lot of time for evaluating a series of let's say 20 observations. I think the slowness of the integration is related to the fact that the p(μ,σ) is calculated using the product of two log-normal distribution but I simply use linear (not exponential) interval for the integrands. Commented Feb 18, 2018 at 16:40