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 dustributiondistribution $N(\mu,\sigma)$ is zero (or I'm wrong?).
In addition, what is the best method to calculate (or estimate) this distribution in python?