Let $(\Omega, \mathcal{A}, \mathbb{P})$ be our respective probability space and $X$ a two-dimensional random variable with values on $[0,1]^2$ and probability density function $f(x_1,x_2)$ for $x_1, x_2 \in [0,1]$.
What are the most efficient methods (from numerics) to approximate the cumulative distribution function \begin{equation} F(x, y) = \int_{[0,y]} \int_{[0,x]} f(x_1, x_2) \ dx_1 dx_2 \end{equation} for any $x,y \in [0,1]$ when a direct computation of the integral is not possible? Note I am interested in evaluating $F(x,y)$ for different values.