I need to find the probability of sampling a specific point on a Gaussian Distribution.
The catch is that the mean of the first Gaussian Distribution is itself sampled from a Gaussian Distribution.
The standard deviation is known, and the mean of the second distribution is 0.
So, calling the variable to sample y:
$ x \sim N(0, \rho) $
$ y\mid x \sim N(x, \sigma) $
To find out the probability of sampling y, I believe I should resolve:
$$ p(y) = \int^{+\infty}_{-\infty} p(y \mid x) \, p(x) \;\mathrm{d}x= \int^{+\infty}_{-\infty} \frac{1}{\sigma \sqrt{2 \pi}} \mathrm{e}^{- \frac{(y-x)^2}{2\sigma^2}} \frac{1}{\rho \sqrt{2 \pi}} \mathrm{e}^{- \frac{x^2}{2\rho^2}}\;\mathrm{d}x $$
which is the probability of sampling y from a Gaussian centered in x times the probability of sampling x from the other Gaussian, for every x.
I tried resolving it on wolframalpha but unfortunately it only works when you specify the standard deviations, while I'd like to have it parametric ( the experession will be evaluated in a program)
I also looked on the internet, but I only found situations in which both the normal distributions depended on the same variable.
I started deriving it by hand, but my calculus exam was long ago and I'm not sure I can do it correctly anymore and unfortunately I don't know any numerical computation program (like Matlab, Scilab, etc)
How can I resolve the given integral to obtain an expression which only depends on y, $\sigma$ and $\rho$?
The value of the double defined integral is : $\quad\frac{1}{\sqrt{2\pi (\sigma^2+\rho^2) }}e^{-\frac{y^2}{2(\sigma^2+\rho^2)}}$