Likelihood for the product of two Gaussian random variables

Questions
1

I have 100 observations $[z_0,…z_99]$ of a random variable $Z=X*Y$ which is a product of two normal random variables X and Y where $X~N(\mu_1,\sigma_1)$ and $Y~N(\mu_2,\sigma_2)$ with contraints $0<X<1$ and $0<Y<1$.

I derived an integral representation for the pdf of Z but there is no analytical formula for the likelihood of Z.

Is there any trick to write likelihood of Z in pymc3 and avoid numerical calculation of cumbersome integrals. Something like:
pm.likelihood_of_z(‘Z’,\mu_1,\sigma_1,\mu_2,\sigma_2,observed= [z_0,…z_99])

with a goal to use MCMC to find $\mu_1,\sigma_1,\mu_2,\sigma_2$.

2

Maybe you can try using this integral directly? eg: Custom theano Op to do numerical integration - #6 by aseyboldt

3

junpenglao, This is a different case as you have an integral representation for the location parameter $\mu$ while the likelihood is still Gaussian.

In my case, i need to write the likelihood similar to formula 22 in http://pldml.icm.edu.pl/pldml/element/bwmeta1.element.bwnjournal-article-doi-10_7151_dmps_1146/c/dmps.1146.pdf
where $mu_1$ and $\mu_2$ are pymc3 random variables. It is not clear how to do it.