@Override
    public double cdf(double[] x) {
        if (x.length != dim) {
            throw new IllegalArgumentException("Sample has different dimension.");
        }

        int Nmax = 10000;
        double alph = GaussianDistribution.getInstance().quantile(0.999);
        double errMax = 0.001;

        double[] v = x.clone();
        MathEx.sub(v, mu);

        double p = 0.0;
        double varSum = 0.0;

        // d is always zero
        double[] e = new double[dim];
        double[] f = new double[dim];
        e[0] = GaussianDistribution.getInstance().cdf(v[0] / sigmaL.get(0, 0));
        f[0] = e[0];

        double[] y = new double[dim];

        double err = 2 * errMax;
        int N;
        for (N = 1; err > errMax && N <= Nmax; N++) {
            double[] w = MathEx.random(dim - 1);
            for (int i = 1; i < dim; i++) {
                y[i - 1] = GaussianDistribution.getInstance().quantile(w[i - 1] * e[i - 1]);
                double q = 0.0;
                for (int j = 0; j < i; j++) {
                    q += sigmaL.get(i, j) * y[j];
                }

                e[i] = GaussianDistribution.getInstance().cdf((v[i] - q) / sigmaL.get(i, i));
                f[i] = e[i] * f[i - 1];
            }

            double del = (f[dim - 1] - p) / N;
            p += del;
            varSum = (N - 2) * varSum / N + del * del;
            err = alph * Math.sqrt(varSum);
        }

        return p;
    }
@Test
    public void testCdf2() {
        System.out.println("cdf2");
        double[][] S = {
                {3.260127902272362, 2.343938296424249, 0.1409050254343716, -0.1628775438743266},
                {2.343938296424249, 4.213034991388330, 1.3997210599608563,  0.3373448510018783},
                {0.1409050254343716, 1.3997210599608563, 4.6042485263677939,  0.0807267064408651},
                {-0.1628775438743266, 0.3373448510018783, 0.0807267064408651,  5.4950949215890672}
        };

        double[] M = {-0.683477474844462,  1.480296478403701,  1.008431991316523,  0.448404211078558};
        double[] X = {0.713919336274493, 0.584408785741822, 0.263119200077829, 0.732513610871908};

        MultivariateGaussianDistribution instance = new MultivariateGaussianDistribution(M, Matrix.of(S));

        // The expected value is based on R
        assertEquals(0.0904191282120575, instance.cdf(X), 1E-3);
    }