private static DoubleUnaryOperator createBinomialModel()

in commons-statistics-inference/src/main/java/org/apache/commons/statistics/inference/UnconditionedExactTest.java [1006:1072]
52 lines of code
9 McCabe index (conditional complexity)

    private static DoubleUnaryOperator createBinomialModel(XYList tableList) {
        final int m = tableList.getMaxX();
        final int n = tableList.getMaxY();
        final int mn = m + n;
        // Compute the probability using logs
        final double[] c = new double[tableList.size()];
        final int[] ij = new int[tableList.size()];
        final int width = tableList.getWidth();

        // Compute the log binomial dynamically for a small number of values
        final IntToDoubleFunction binomM;
        final IntToDoubleFunction binomN;
        if (tableList.size() < mn) {
            binomM = k -> LogBinomialCoefficient.value(m, k);
            binomN = k -> LogBinomialCoefficient.value(n, k);
        } else {
            // Pre-compute all values
            binomM = createLogBinomialCoefficients(m);
            binomN = m == n ? binomM : createLogBinomialCoefficients(n);
        }

        // Handle special cases i+j == 0 and i+j == m+n.
        // These will occur only once, if at all. Mark if they occur.
        int flag = 0;
        int j = 0;
        for (int i = 0; i < c.length; i++) {
            final int index = tableList.get(i);
            final int x = index % width;
            final int y = index / width;
            final int xy = x + y;
            if (xy == 0) {
                flag |= 1;
            } else if (xy == mn) {
                flag |= 2;
            } else {
                ij[j] = xy;
                c[j] = binomM.applyAsDouble(x) + binomN.applyAsDouble(y);
                j++;
            }
        }

        final int size = j;
        final boolean ij0 = (flag & 1) != 0;
        final boolean ijmn = (flag & 2) != 0;
        return pi -> {
            final double logp = Math.log(pi);
            final double log1mp = Math.log1p(-pi);
            double sum = 0;
            for (int i = 0; i < size; i++) {
                // binom(m, i) * binom(n, j) * pi^(i+j) * (1-pi)^(m+n-i-j)
                sum += Math.exp(ij[i] * logp + (mn - ij[i]) * log1mp + c[i]);
            }
            // Add the simplified terms where the binomial is 1.0 and one power is x^0 == 1.0.
            // This avoids 0 * log(x) generating NaN when x is 0 in the case where pi was 0 or 1.
            // Reuse exp (not pow) to support pi approaching 0 or 1.
            if (ij0) {
                // pow(1-pi, mn)
                sum += Math.exp(mn * log1mp);
            }
            if (ijmn) {
                // pow(pi, mn)
                sum += Math.exp(mn * logp);
            }
            // The optimizer minimises the function so this returns -p.
            return -sum;
        };
    }