private static double moment2()

in commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/TruncatedNormalDistribution.java [469:545]
30 lines of code
5 McCabe index (conditional complexity)

    private static double moment2(double a, double b) {
        // Assume a < b.
        // a == b is handled in the variance method
        if (Math.abs(a) > Math.abs(b)) {
            return moment2(-b, -a);
        }

        // Here:
        // |a| <= |b|
        // a < b
        // 0 < b

        if (a <= -MAX_X) {
            // No truncation
            return 1;
        }
        if (b >= MAX_X) {
            // One-sided truncation.
            // For a -> inf : moment2 -> a*a
            // This occurs when erfcx(z) is approximated by (1/sqrt(pi)) / z and terms
            // cancel. z > 6.71e7, a > 9.49e7
            return 1 + ROOT_2_PI * a / Erfcx.value(a / ROOT2);
        }

        // pdf = exp(-0.5*x*x) / sqrt(2*pi)
        // cdf = erfc(-x/sqrt(2)) / 2
        // Compute:
        // 1 - (b*pdf(b) - a*pdf(a)) / cdf(b, a)
        // = (cdf(b, a) - b*pdf(b) -a*pdf(a)) / cdf(b, a)

        // Note:
        // For z -> 0:
        //   sqrt(pi / 2) * erf(z / sqrt(2)) -> z
        //   z * Math.exp(-0.5 * z * z) -> z
        // Both computations below have cancellation as b -> 0 and the
        // second moment is not computable as the fraction P/Q
        // since P < ulp(Q). This always occurs when b < MIN_X
        // if MIN_X is set at the point where
        //   exp(-0.5 * z * z) / sqrt(2 pi) == 1 / sqrt(2 pi).
        // This is JDK dependent due to variations in Math.exp.
        // For b < MIN_X the second moment can be approximated using
        // a uniform distribution: (b^3 - a^3) / (3b - 3a).
        // In practice it also occurs when b > MIN_X since any a < MIN_X
        // is effectively zero for part of the computation. A
        // threshold to transition to a uniform distribution
        // approximation is a compromise. Also note it will not
        // correct computation when (b-a) is small and is far from 0.
        // Thus the second moment is left to be inaccurate for
        // small ranges (b-a) and the variance -> 0 when the true
        // variance is close to or below machine epsilon.

        double m;

        if (a <= 0) {
            // Opposite signs
            final double ea = ROOT_PI_2 * Erf.value(a / ROOT2);
            final double eb = ROOT_PI_2 * Erf.value(b / ROOT2);
            final double fa = ea - a * Math.exp(-0.5 * a * a);
            final double fb = eb - b * Math.exp(-0.5 * b * b);
            // Assume fb >= fa && eb >= ea
            // If fb <= fa this is a tiny range around 0
            m = (fb - fa) / (eb - ea);
            // Clip to the range
            m = clip(m, 0, 1);
        } else {
            final double dx = 0.5 * (b + a) * (b - a);
            final double ex = Math.exp(-dx);
            final double ea = ROOT_PI_2 * Erfcx.value(a / ROOT2);
            final double eb = ROOT_PI_2 * Erfcx.value(b / ROOT2);
            final double fa = ea + a;
            final double fb = eb + b;
            m = (fa - fb * ex) / (ea - eb * ex);
            // Clip to the range
            m = clip(m, a * a, b * b);
        }
        return m;
    }