integration – Deterministic Approximation of Products of Integrals

Suppose that I am interested in approximating the product of two integrals. For simplicity, let $P$ be a probability density function, let $f, g$ be two $P$-integrable functions, define

phi = int P(x) f(x) , dx \
gamma = int P(x) g(x) , dx,

and say that I am interested in approximating the product $phi cdot gamma$.

I am interested in how a numerical analyst would approach this problem. A naive solution to this would be to do e.g. Gaussian quadrature with $P$ as the base measure, use the GQ nodes to form approximations of $phi, gamma$, and then multiply the approximations together. Can / does one ever do fancier things than this in practice?

I ask because I come from a statistics / applied probability background, where one often places a premium on having unbiased estimates of integrals. As such, taking the product of two estimates which use the same set of function evaluations will generally induce a statistical bias, and so there are simple and popular techniques which allow for this bias to be removed. I am curious about whether similar approaches are considered in numerical analysis.

To compress the question further: from the point of view of numerical analysis, are there better ways of approximating a product of integrals than taking the product of approximations to each of the integrals?