I would like to know if there exists a way in which the double integral
$$\int_0^{\infty}d\tau f(\tau)\int_0^{\tau}d\theta g(\tau-\theta)h(\theta)$$
can be computed (numerically) efficiently. It is similar to the integral of the convolution between $g$ and $h$ weighted by the function $f$, and I was wondering if I could exploit a result akin to this one.
Some additional details:
In reality the function $h(t)$ is $h(x,t)$, and it's a p.d.f over $x$ for every value of $t$.
$$\int_{-\infty}^{\infty} dx\, h(x,t) = 1 \quad \forall t$$
Also the functions $f$ and $g$ are p.d.f.s so they satisfy
$$\int_0^{\infty} dt\, f(t) = \int_0^{\infty} dt\, g(t) = 1$$
The domain of $f(t)$, $g(t)$ and $h(\cdot,t)$ is $t \in [0,\infty)$