Numerical CDF estimation for complicated random variable Given a combination $U$ of several random variables $X,Y,Z...$ with known distributions, what is an efficient numerical algorithm to estimate PDF or CDF of $U$, if its CDF, PDF, characteristic function cannot be expressed analytically, and at the points of interest the PDF is very small? The direct Monte Carlo sampling would be very inaccurate at those points.
As a simple example, how to estimate the CDF $F(t)=\Pr[U<t]$ of $U=(X+Y)(X+Z)$, where $X,Y,Z$ are iid standard Gaussians, at $t=-10$?
 A: The best solution would be to derive the exact form of the pdf of $U = (X+Y)(X+Z)$. But if the question is more general, or if closed-form solutions are not available, then the way I would proceed would be to:
Step 1:  Define the joint pdf 
Define the joint pdf of $(X,Y,Z)$, say $f(x,y,z)$ ... which is simply the product of the 3 standard Normal pdfs:

      f = Exp[-x^2/2 - y^2/2 - z^2/2]/(2 Sqrt[2] Pi^(3/2));

domain[f] = {{x, -Infinity, Infinity}, {y, -Infinity, Infinity}, {z, -Infinity, Infinity}};

Step 2:  numerically integrate
To find any arbitrary probability, numerically integrate.
In your case, we seek: $P\big((X+Y)(X+Z) < - 10\big)$.
I am using Mathematica syntax here, but any decent package that supports arbitrary precision numerics should do the job:
NIntegrate[ Boole[(x+y)(x+z) < -10] * f ,
{x, -Infinity, Infinity}, {y, -Infinity, Infinity}, {z, -Infinity, Infinity}]
which returns:

3.80999*10^-6

One can impose desired working precisions and so on.  The Boole function operates like an indicator function:  this allows one to calculate the probability of essentially any desired expression, so the approach is very flexible.
