static double moment1()

in commons-statistics-distribution/src/main/java/org/apache/commons/statistics/distribution/TruncatedNormalDistribution.java [407:458]

• 26 lines of code
• 7 McCabe index (conditional complexity)

    static double moment1(double a, double b) {
        // Assume a <= b
        if (a == b) {
            return a;
        }
        if (Math.abs(a) > Math.abs(b)) {
            // Subtract from zero to avoid generating -0.0
            return 0 - moment1(-b, -a);
        }

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

        if (a <= -MAX_X) {
            // No truncation
            return 0;
        }
        if (b >= MAX_X) {
            // One-sided truncation
            return ROOT_2_PI / Erfcx.value(a / ROOT2);
        }

        // pdf = exp(-0.5*x*x) / sqrt(2*pi)
        // cdf = erfc(-x/sqrt(2)) / 2
        // Compute:
        // -(pdf(b) - pdf(a)) / cdf(b, a)
        // Note:
        // exp(-0.5*b*b) - exp(-0.5*a*a)
        // Use cancellation of powers:
        // exp(-0.5*(b*b-a*a)) * exp(-0.5*a*a) - exp(-0.5*a*a)
        // expm1(-0.5*(b*b-a*a)) * exp(-0.5*a*a)

        // dx = -0.5*(b*b-a*a)
        final double dx = 0.5 * (b + a) * (b - a);
        final double m;
        if (a <= 0) {
            // Opposite signs
            m = ROOT_2_PI * -Math.expm1(-dx) * Math.exp(-0.5 * a * a) / ErfDifference.value(a / ROOT2, b / ROOT2);
        } else {
            final double z = Math.exp(-dx) * Erfcx.value(b / ROOT2) - Erfcx.value(a / ROOT2);
            if (z == 0) {
                // Occurs when a and b have large magnitudes and are very close
                return (a + b) * 0.5;
            }
            m = ROOT_2_PI * Math.expm1(-dx) / z;
        }

        // Clip to the range
        return clip(m, a, b);
    }