private static SharedStateDiscreteSampler createBinomialDistributionSampler()

in commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/MarsagliaTsangWangDiscreteSampler.java [1049:1095]


        private static SharedStateDiscreteSampler createBinomialDistributionSampler(
                UniformRandomProvider rng, int trials, double probabilityOfSuccess) {

            // The maximum supported value for Math.exp is approximately -744.
            // This occurs when trials is large and p is close to 1.
            // Handle this by using an inversion: generate j=Binomial(n,1-p), return n-j
            final boolean useInversion = probabilityOfSuccess > 0.5;
            final double p = useInversion ? 1 - probabilityOfSuccess : probabilityOfSuccess;

            // Check if the distribution can be computed
            final double p0 = Math.exp(trials * Math.log(1 - p));
            if (p0 < Double.MIN_VALUE) {
                throw new IllegalArgumentException("Unable to compute distribution");
            }

            // First find size of probability array
            double t = p0;
            final double h = p / (1 - p);
            // Find first probability above the threshold of 2^-31
            int begin = 0;
            if (t * DOUBLE_31 < 1) {
                // Somewhere after p(0)
                // Note:
                // If this loop is entered p(0) is < 2^-31.
                // This has been tested at the extreme for p(0)=Double.MIN_VALUE and either
                // p=0.5 or trials=2^16-1 and does not fail to find the beginning.
                for (int i = 1; i <= trials; i++) {
                    t *= (trials + 1 - i) * h / i;
                    if (t * DOUBLE_31 >= 1) {
                        begin = i;
                        break;
                    }
                }
            }
            // Find last probability
            int end = trials;
            for (int i = begin + 1; i <= trials; i++) {
                t *= (trials + 1 - i) * h / i;
                if (t * DOUBLE_31 < 1) {
                    end = i - 1;
                    break;
                }
            }

            return createBinomialDistributionSamplerFromRange(rng, trials, p, useInversion,
                    p0, begin, end);
        }