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:
with a goal to use MCMC to find $\mu_1,\sigma_1,\mu_2,\sigma_2$.