/*
 * Licensed to the Apache Software Foundation (ASF) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * The ASF licenses this file to You under the Apache License, Version 2.0
 * (the "License"); you may not use this file except in compliance with
 * the License.  You may obtain a copy of the License at
 *
 *      http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
package org.apache.commons.statistics.distribution;

import java.util.Arrays;

import org.apache.commons.rng.simple.RandomSource;
import org.apache.commons.rng.UniformRandomProvider;
import org.apache.commons.math3.stat.StatUtils;
import org.apache.commons.math3.stat.inference.GTest;
import org.junit.jupiter.api.Assertions;
import org.junit.jupiter.api.Test;

class BetaDistributionTest {

    static final double[] ALPHA_BETAS = {0.1, 1, 10, 100, 1000};
    static final double EPSILON = StatUtils.min(ALPHA_BETAS);
    static final double[] CUMULATIVE_TEST_POINTS = new double[]{
        -0.1, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1
    };

    @Test
    void testCumulative() {
        final double[] x = CUMULATIVE_TEST_POINTS;
        // all test data computed using R 2.5
        checkCumulative(0.1, 0.1,
                        x, new double[]{
                            0.0000000000, 0.0000000000, 0.4063850939, 0.4397091902, 0.4628041861,
                            0.4821200456, 0.5000000000, 0.5178799544, 0.5371958139, 0.5602908098,
                            0.5936149061, 1.0000000000, 1.0000000000});
        checkCumulative(0.1, 0.5,
                        x, new double[]{
                            0.0000000000, 0.0000000000, 0.7048336221, 0.7593042194, 0.7951765304,
                            0.8234948385, 0.8480017124, 0.8706034370, 0.8926585878, 0.9156406404,
                            0.9423662883, 1.0000000000, 1.0000000000});
        checkCumulative(0.1, 1.0,
                        x, new double[]{
                            0.0000000000, 0.0000000000, 0.7943282347, 0.8513399225, 0.8865681506,
                            0.9124435366, 0.9330329915, 0.9502002165, 0.9649610951, 0.9779327685,
                            0.9895192582, 1.0000000000, 1.0000000000});
        checkCumulative(0.1, 2.0,
                        x, new double[]{
                            0.0000000000, 0.0000000000, 0.8658177758, 0.9194471163, 0.9486279211,
                            0.9671901487, 0.9796846411, 0.9882082252, 0.9939099280, 0.9974914239,
                            0.9994144508, 1.0000000000, 1.0000000000});
        checkCumulative(0.1, 4.0,
                        x, new double[]{
                            0.0000000000, 0.0000000000, 0.9234991121, 0.9661958941, 0.9842285085,
                            0.9928444112, 0.9970040660, 0.9989112804, 0.9996895625, 0.9999440793,
                            0.9999967829, 1.0000000000, 1.0000000000});
        checkCumulative(0.5, 0.1,
                        x, new double[]{
                            0.00000000000, 0.00000000000, 0.05763371168, 0.08435935962,
                            0.10734141216, 0.12939656302, 0.15199828760, 0.17650516146,
                            0.20482346963, 0.24069578055, 0.29516637795, 1.00000000000, 1.00000000000});

        checkCumulative(0.5, 0.5,
                        x, new double[]{
                            0.0000000000, 0.0000000000, 0.2048327647, 0.2951672353, 0.3690101196,
                            0.4359057832, 0.5000000000, 0.5640942168, 0.6309898804, 0.7048327647,
                            0.7951672353, 1.0000000000, 1.0000000000});
        checkCumulative(0.5, 1.0,
                        x, new double[]{
                            0.0000000000, 0.0000000000, 0.3162277660, 0.4472135955, 0.5477225575,
                            0.6324555320, 0.7071067812, 0.7745966692, 0.8366600265, 0.8944271910,
                            0.9486832981, 1.0000000000, 1.0000000000});
        checkCumulative(0.5, 2.0,
                        x, new double[]{
                            0.0000000000, 0.0000000000, 0.4585302607, 0.6260990337, 0.7394254526,
                            0.8221921916, 0.8838834765, 0.9295160031, 0.9621590305, 0.9838699101,
                            0.9961174630, 1.0000000000, 1.0000000000});
        checkCumulative(0.5, 4.0,
                        x, new double[]{
                            0.0000000000, 0.0000000000, 0.6266250826, 0.8049844719, 0.8987784842,
                            0.9502644369, 0.9777960959, 0.9914837366, 0.9974556254, 0.9995223859,
                            0.9999714889, 1.0000000000, 1.0000000000});
        checkCumulative(1.0, 0.1,
                        x, new double[]{
                            0.00000000000, 0.00000000000, 0.01048074179, 0.02206723146,
                            0.03503890488, 0.04979978349, 0.06696700846, 0.08755646344,
                            0.11343184943, 0.14866007748, 0.20567176528, 1.00000000000, 1.00000000000});
        checkCumulative(1.0, 0.5,
                        x, new double[]{
                            0.00000000000, 0.00000000000, 0.05131670195, 0.10557280900,
                            0.16333997347, 0.22540333076, 0.29289321881, 0.36754446797,
                            0.45227744249, 0.55278640450, 0.68377223398, 1.00000000000, 1.00000000000});
        checkCumulative(1, 1,
                        x, new double[]{
                            0.0, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.0});
        checkCumulative(1, 2,
                        x, new double[]{
                            0.00, 0.00, 0.19, 0.36, 0.51, 0.64, 0.75, 0.84, 0.91, 0.96, 0.99, 1.00, 1.00});
        checkCumulative(1, 4,
                        x, new double[]{
                            0.0000, 0.0000, 0.3439, 0.5904, 0.7599, 0.8704, 0.9375, 0.9744, 0.9919,
                            0.9984, 0.9999, 1.0000, 1.0000});
        checkCumulative(2.0, 0.1,
                        x, new double[]{
                            0.0000000000000, 0.0000000000000, 0.0005855492117, 0.0025085760862,
                            0.0060900720266, 0.0117917748341, 0.0203153588864, 0.0328098512512,
                            0.0513720788952, 0.0805528836776, 0.1341822241505, 1.0000000000000, 1.0000000000000});
        checkCumulative(2, 1,
                        x, new double[]{
                            0.00, 0.00, 0.01, 0.04, 0.09, 0.16, 0.25, 0.36, 0.49, 0.64, 0.81, 1.00, 1.00});
        checkCumulative(2.0, 0.5,
                        x, new double[]{
                            0.000000000000, 0.000000000000, 0.003882537047, 0.016130089900,
                            0.037840969486, 0.070483996910, 0.116116523517, 0.177807808356,
                            0.260574547368, 0.373900966300, 0.541469739276, 1.000000000000, 1.000000000000});
        checkCumulative(2, 2,
                        x, new double[]{
                            0.000, 0.000, 0.028, 0.104, 0.216, 0.352, 0.500, 0.648, 0.784, 0.896, 0.972, 1.000, 1.000});
        checkCumulative(2, 4,
                        x, new double[]{
                            0.00000, 0.00000, 0.08146, 0.26272, 0.47178, 0.66304, 0.81250, 0.91296,
                            0.96922, 0.99328, 0.99954, 1.00000, 1.00000});
        checkCumulative(4.0, 0.1,
                        x, new double[]{
                            0.000000000e+00, 0.000000000e+00, 3.217128269e-06, 5.592070271e-05,
                            3.104375474e-04, 1.088719595e-03, 2.995933981e-03, 7.155588777e-03,
                            1.577149153e-02, 3.380410585e-02, 7.650088789e-02, 1.000000000e+00, 1.000000000e+00});
        checkCumulative(4.0, 0.5,
                        x, new double[]{
                            0.000000000e+00, 0.000000000e+00, 2.851114863e-05, 4.776140576e-04,
                            2.544374616e-03, 8.516263371e-03, 2.220390414e-02, 4.973556312e-02,
                            1.012215158e-01, 1.950155281e-01, 3.733749174e-01, 1.000000000e+00, 1.000000000e+00});
        checkCumulative(4, 1,
                        x, new double[]{
                            0.0000, 0.0000, 0.0001, 0.0016, 0.0081, 0.0256, 0.0625, 0.1296, 0.2401,
                            0.4096, 0.6561, 1.0000, 1.0000});
        checkCumulative(4, 2,
                        x, new double[]{
                            0.00000, 0.00000, 0.00046, 0.00672, 0.03078, 0.08704, 0.18750, 0.33696,
                            0.52822, 0.73728, 0.91854, 1.00000, 1.00000});
        checkCumulative(4, 4,
                        x, new double[]{
                            0.000000, 0.000000, 0.002728, 0.033344, 0.126036, 0.289792, 0.500000,
                            0.710208, 0.873964, 0.966656, 0.997272, 1.000000, 1.000000});
    }

    @Test
    void testSurvival() {
        checkSurvival(1, 1, new double[]{
            0.0, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.0});
        checkSurvival(4, 1, new double[]{
            0.0000, 0.0000, 0.0001, 0.0016, 0.0081, 0.0256, 0.0625, 0.1296, 0.2401,
            0.4096, 0.6561, 1.0000, 1.0000});
        checkSurvival(2, 4, new double[]{
            0.00000, 0.00000, 0.08146, 0.26272, 0.47178, 0.66304, 0.81250, 0.91296,
            0.96922, 0.99328, 0.99954, 1.00000, 1.00000});
    }

    @Test
    void testCumulativeAndSurvivalComplement() {
        checkCumulativeSurvivalComplement(0.1, 0.1);
        checkCumulativeSurvivalComplement(0.1, 0.5);
        checkCumulativeSurvivalComplement(0.1, 1.0);
        checkCumulativeSurvivalComplement(0.1, 2.0);
        checkCumulativeSurvivalComplement(0.1, 4.0);
        checkCumulativeSurvivalComplement(0.5, 0.1);
        checkCumulativeSurvivalComplement(0.5, 0.5);
        checkCumulativeSurvivalComplement(0.5, 1.0);
        checkCumulativeSurvivalComplement(0.5, 2.0);
        checkCumulativeSurvivalComplement(0.5, 4.0);
        checkCumulativeSurvivalComplement(1.0, 0.1);
        checkCumulativeSurvivalComplement(1.0, 0.5);
        checkCumulativeSurvivalComplement(1, 1);
        checkCumulativeSurvivalComplement(1, 2);
        checkCumulativeSurvivalComplement(1, 4);
        checkCumulativeSurvivalComplement(2.0, 0.1);
        checkCumulativeSurvivalComplement(2, 1);
        checkCumulativeSurvivalComplement(2.0, 0.5);
        checkCumulativeSurvivalComplement(2, 2);
        checkCumulativeSurvivalComplement(2, 4);
        checkCumulativeSurvivalComplement(4.0, 0.1);
        checkCumulativeSurvivalComplement(4.0, 0.5);
        checkCumulativeSurvivalComplement(4, 1);
        checkCumulativeSurvivalComplement(4, 2);
        checkCumulativeSurvivalComplement(4, 4);
    }

    /** Precision tests for verifying that CDF calculates accurately in cases where 1-cdf(x) is inaccurately 1. */
    @Test
    void testCumulativePrecision() {
        // Calculated using WolframAlpha
        checkCumulativePrecision(5.0, 5.0, 0.0001, 1.2595800539968654e-18);
        checkCumulativePrecision(4.0, 5.0, 0.00001, 6.999776002800025e-19);
        checkCumulativePrecision(5.0, 4.0, 0.0001, 5.598600119996539e-19);
        checkCumulativePrecision(6.0, 2.0, 0.001, 6.994000000000028e-18);
        checkCumulativePrecision(2.0, 6.0, 1e-9, 2.0999999930000014e-17);
    }

    /**
     * Precision tests for verifying that survivalFunction calculates accurately in cases
     * where 1-sf(x) is inaccurately 1.
     */
    @Test
    void testSurvivalPrecision() {
        // Calculated using WolframAlpha
        checkSurvivalPrecision(5.0, 5.0, 0.9999, 1.2595800539961496e-18);
        checkSurvivalPrecision(4.0, 5.0, 0.9999, 5.598600119993397e-19);
        checkSurvivalPrecision(5.0, 4.0, 0.99998, 1.1199283217964632e-17);
        checkSurvivalPrecision(6.0, 2.0, 0.999999999, 2.0999998742158932e-17);
        checkSurvivalPrecision(2.0, 6.0, 0.999, 6.994000000000077e-18);
    }

    private static void checkCumulative(double alpha, double beta, double[] x, double[] cumes) {
        final BetaDistribution dist = new BetaDistribution(alpha, beta);
        for (int i = 0; i < x.length; i++) {
            Assertions.assertEquals(cumes[i], dist.cumulativeProbability(x[i]), 1e-8);
        }

        for (int i = 1; i < x.length - 1; i++) {
            Assertions.assertEquals(x[i], dist.inverseCumulativeProbability(cumes[i]), 1e-5);
        }
    }

    private static void checkSurvival(double alpha, double beta, double[] cumes) {
        final BetaDistribution dist = new BetaDistribution(alpha, beta);
        for (int i = 0; i < CUMULATIVE_TEST_POINTS.length; i++) {
            Assertions.assertEquals(1 - cumes[i], dist.survivalProbability(CUMULATIVE_TEST_POINTS[i]), 1e-8);
        }
    }

    private static void checkCumulativeSurvivalComplement(double alpha, double beta) {
        final BetaDistribution dist = new BetaDistribution(alpha, beta);
        for (int i = 0; i < CUMULATIVE_TEST_POINTS.length; i++) {
            final double x = CUMULATIVE_TEST_POINTS[i];
            Assertions.assertEquals(1, dist.cumulativeProbability(x) + dist.survivalProbability(x), 1e-8);
        }
    }

    private static void checkCumulativePrecision(double alpha, double beta, double value, double expected) {
        final double tolerance = 1e-22;
        final BetaDistribution dist = new BetaDistribution(alpha, beta);
        Assertions.assertEquals(
            dist.cumulativeProbability(value),
            expected,
            tolerance,
            () -> "cumulative probability not precise at " + value + " for a=" + alpha + " & b=" + beta);
    }

    private static void checkSurvivalPrecision(double alpha, double beta, double value, double expected) {
        final double tolerance = 1e-22;
        final BetaDistribution dist = new BetaDistribution(alpha, beta);
        Assertions.assertEquals(
            dist.survivalProbability(value),
            expected,
            tolerance,
            () -> "survival function not precise at " + value + " for a=" + alpha + " & b=" + beta);
    }

    @Test
    void testDensity() {
        final double[] x = new double[]{1e-6, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9};
        checkDensity(0.1, 0.1,
                     x, new double[]{
                         12741.2357380649, 0.4429889586665234, 2.639378715e-01, 2.066393611e-01,
                         1.832401831e-01, 1.766302780e-01, 1.832404579e-01, 2.066400696e-01,
                         2.639396531e-01, 4.429925026e-01});
        checkDensity(0.1, 0.5,
                     x, new double[]{
                         2.218377102e+04, 7.394524202e-01, 4.203020268e-01, 3.119435533e-01,
                         2.600787829e-01, 2.330648626e-01, 2.211408259e-01, 2.222728708e-01,
                         2.414013907e-01, 3.070567405e-01});
        checkDensity(0.1, 1.0,
                     x, new double[]{
                         2.511886432e+04, 7.943210858e-01, 4.256680458e-01, 2.955218303e-01,
                         2.281103709e-01, 1.866062624e-01, 1.583664652e-01, 1.378514078e-01,
                         1.222414585e-01, 1.099464743e-01});
        checkDensity(0.1, 2.0,
                     x, new double[]{
                         2.763072312e+04, 7.863770012e-01, 3.745874120e-01, 2.275514842e-01,
                         1.505525939e-01, 1.026332391e-01, 6.968107049e-02, 4.549081293e-02,
                         2.689298641e-02, 1.209399123e-02});
        checkDensity(0.1, 4.0,
                     x, new double[]{
                         2.997927462e+04, 6.911058917e-01, 2.601128486e-01, 1.209774010e-01,
                         5.880564714e-02, 2.783915474e-02, 1.209657335e-02, 4.442148268e-03,
                         1.167143939e-03, 1.312171805e-04});
        checkDensity(0.5, 0.1,
                     x, new double[]{
                         88.3152184726, 0.3070542841, 0.2414007269, 0.2222727015,
                         0.2211409364, 0.2330652355, 0.2600795198, 0.3119449793,
                         0.4203052841, 0.7394649088});
        checkDensity(0.5, 0.5,
                     x, new double[]{
                         318.3100453389, 1.0610282383, 0.7957732234, 0.6946084565,
                         0.6497470636, 0.6366197724, 0.6497476051, 0.6946097796,
                         0.7957762075, 1.0610376697});
        checkDensity(0.5, 1.0,
                     x, new double[]{
                         500.0000000000, 1.5811309244, 1.1180311937, 0.9128694077,
                         0.7905684268, 0.7071060741, 0.6454966865, 0.5976138778,
                         0.5590166450, 0.5270459839});
        checkDensity(0.5, 2.0,
                     x, new double[]{
                         749.99925000000, 2.134537420613655, 1.34163575536, 0.95851150881,
                         0.71151039830, 0.53032849490, 0.38729704363, 0.26892534859,
                         0.16770415497, 0.07905610701});
        checkDensity(0.5, 4.0,
                     x, new double[]{
                         1.093746719e+03, 2.52142232809988, 1.252190241e+00, 6.849343920e-01,
                         3.735417140e-01, 1.933481570e-01, 9.036885833e-02, 3.529621669e-02,
                         9.782644546e-03, 1.152878503e-03});
        checkDensity(1.0, 0.1,
                     x, new double[]{
                         0.1000000900, 0.1099466942, 0.1222417336, 0.1378517623, 0.1583669403,
                         0.1866069342, 0.2281113974, 0.2955236034, 0.4256718768,
                         0.7943353837});
        checkDensity(1.0, 0.5,
                     x, new double[]{
                         0.5000002500, 0.5270465695, 0.5590173438, 0.5976147315, 0.6454977623,
                         0.7071074883, 0.7905704033, 0.9128724506,
                         1.1180367838, 1.5811467358});
        checkDensity(1, 1,
                     x, new double[]{
                         1, 1, 1,
                         1, 1, 1, 1,
                         1, 1, 1});
        checkDensity(1, 2,
                     x, new double[]{
                         1.999998, 1.799998, 1.599998, 1.399998, 1.199998, 0.999998, 0.799998,
                         0.599998, 0.399998,
                         0.199998});
        checkDensity(1, 4,
                     x, new double[]{
                         3.999988000012, 2.915990280011, 2.047992320010, 1.371994120008,
                         0.863995680007, 0.499997000006, 0.255998080005, 0.107998920004,
                         0.031999520002, 0.003999880001});
        checkDensity(2.0, 0.1,
                     x, new double[]{
                         1.100000990e-07, 1.209425730e-02, 2.689331586e-02, 4.549123318e-02,
                         6.968162794e-02, 1.026340191e-01, 1.505537732e-01, 2.275534997e-01,
                         3.745917198e-01, 7.863929037e-01});
        checkDensity(2.0, 0.5,
                     x, new double[]{
                         7.500003750e-07, 7.905777599e-02, 1.677060417e-01, 2.689275256e-01,
                         3.872996256e-01, 5.303316769e-01, 7.115145488e-01, 9.585174425e-01,
                         1.341645818e+00, 2.134537420613655});
        checkDensity(2, 1,
                     x, new double[]{
                         0.000002, 0.200002, 0.400002, 0.600002, 0.800002, 1.000002, 1.200002,
                         1.400002, 1.600002,
                         1.800002});
        checkDensity(2, 2,
                     x, new double[]{
                         5.9999940e-06, 5.4000480e-01, 9.6000360e-01, 1.2600024e+00,
                         1.4400012e+00, 1.5000000e+00, 1.4399988e+00, 1.2599976e+00,
                         9.5999640e-01, 5.3999520e-01});
        checkDensity(2, 4,
                     x, new double[]{
                         0.00001999994, 1.45800971996, 2.04800255997, 2.05799803998,
                         1.72799567999, 1.24999500000, 0.76799552000, 0.37799676001,
                         0.12799824001, 0.01799948000});
        checkDensity(4.0, 0.1,
                     x, new double[]{
                         1.193501074e-19, 1.312253162e-04, 1.167181580e-03, 4.442248535e-03,
                         1.209679109e-02, 2.783958903e-02, 5.880649983e-02, 1.209791638e-01,
                         2.601171405e-01, 6.911229392e-01});
        checkDensity(4.0, 0.5,
                     x, new double[]{
                         1.093750547e-18, 1.152948959e-03, 9.782950259e-03, 3.529697305e-02,
                         9.037036449e-02, 1.933508639e-01, 3.735463833e-01, 6.849425461e-01,
                         1.252205894e+00, 2.52142232809988});
        checkDensity(4, 1,
                     x, new double[]{
                         4.000000000e-18, 4.000120001e-03, 3.200048000e-02, 1.080010800e-01,
                         2.560019200e-01, 5.000030000e-01, 8.640043200e-01, 1.372005880e+00,
                         2.048007680e+00, 2.916009720e+00});
        checkDensity(4, 2,
                     x, new double[]{
                         1.999998000e-17, 1.800052000e-02, 1.280017600e-01, 3.780032400e-01,
                         7.680044800e-01, 1.250005000e+00, 1.728004320e+00, 2.058001960e+00,
                         2.047997440e+00, 1.457990280e+00});
        checkDensity(4, 4,
                     x, new double[]{
                         1.399995800e-16, 1.020627216e-01, 5.734464512e-01, 1.296547409e+00,
                         1.935364838e+00, 2.187500000e+00, 1.935355162e+00, 1.296532591e+00,
                         5.734335488e-01, 1.020572784e-01});
    }

    @SuppressWarnings("boxing")
    private static void checkDensity(double alpha, double beta, double[] x, double[] expected) {
        final BetaDistribution dist = new BetaDistribution(alpha, beta);
        for (int i = 0; i < x.length; i++) {
            final int index = i;
            Assertions.assertEquals(expected[i], dist.density(x[i]), 1e-5,
                () -> String.format("density at x=%.1f for alpha=%.1f, beta=%.1f", x[index], alpha, beta));
        }
    }

    @Test
    void testDensityExtremes() {
        final BetaDistribution dist = new BetaDistribution(2, 3);
        Assertions.assertEquals(Double.NEGATIVE_INFINITY, dist.logDensity(-0.1));
        Assertions.assertEquals(Double.NEGATIVE_INFINITY, dist.logDensity(1.1));
        Assertions.assertEquals(0.0, dist.density(-0.1));
        Assertions.assertEquals(0.0, dist.density(1.1));

        Assertions.assertEquals(Double.NEGATIVE_INFINITY, dist.logDensity(0.0));
        Assertions.assertEquals(Double.NEGATIVE_INFINITY, dist.logDensity(1.0));
    }

    @Test
    void testLogDensityPrecondition1() {
        final BetaDistribution dist = new BetaDistribution(0.5, 3);
        Assertions.assertThrows(DistributionException.class, () -> dist.logDensity(0.0));
    }

    @Test
    void testLogDensityPrecondition2() {
        final BetaDistribution dist = new BetaDistribution(2, 0.5);
        Assertions.assertThrows(DistributionException.class, () -> dist.logDensity(1.0));
    }

    @Test
    void testMoments() {
        final double tol = 1e-9;
        BetaDistribution dist;

        dist = new BetaDistribution(1, 1);
        Assertions.assertEquals(0.5, dist.getMean(), tol);
        Assertions.assertEquals(1.0 / 12.0, dist.getVariance(), tol);

        dist = new BetaDistribution(2, 5);
        Assertions.assertEquals(2.0 / 7.0, dist.getMean(), tol);
        Assertions.assertEquals(10.0 / (49.0 * 8.0), dist.getVariance(), tol);
    }

    @Test
    void testMomentsSampling() {
        final UniformRandomProvider rng = RandomSource.WELL_1024_A.create(123456789L);
        final int numSamples = 1000;
        for (final double alpha : ALPHA_BETAS) {
            for (final double beta : ALPHA_BETAS) {
                final BetaDistribution betaDistribution = new BetaDistribution(alpha, beta);
                final double[] observed = TestUtils.sample(numSamples,
                        betaDistribution.createSampler(rng));
                Arrays.sort(observed);

                Assertions.assertEquals(betaDistribution.getMean(), StatUtils.mean(observed), EPSILON,
                    () -> String.format("E[Beta(%.2f, %.2f)]", alpha, beta));
                Assertions.assertEquals(betaDistribution.getVariance(), StatUtils.variance(observed), EPSILON,
                    () -> String.format("Var[Beta(%.2f, %.2f)]", alpha, beta));
            }
        }
    }

    @Test
    void testGoodnessOfFit() {
        final UniformRandomProvider rng = RandomSource.WELL_19937_A.create(123456789L);

        final int numSamples = 1000;
        final double level = 0.01;
        for (final double alpha : ALPHA_BETAS) {
            for (final double beta : ALPHA_BETAS) {
                final BetaDistribution betaDistribution = new BetaDistribution(alpha, beta);

                final ContinuousDistribution.Sampler sampler = betaDistribution.createSampler(rng);
                final double[] observed = TestUtils.sample(numSamples, sampler);

                final double gT = gTest(betaDistribution, observed);
                Assertions.assertFalse(gT < level,
                    () -> "G goodness-of-fit (" + gT + ") test rejected null at alpha = " + level);
            }
        }
    }

    private static double gTest(final ContinuousDistribution expectedDistribution, final double[] values) {
        final int numBins = values.length / 30;
        final double[] breaks = new double[numBins];
        for (int b = 0; b < breaks.length; b++) {
            breaks[b] = expectedDistribution.inverseCumulativeProbability((double) b / numBins);
        }

        final long[] observed = new long[numBins];
        for (final double value : values) {
            int b = 0;
            do {
                b++;
            } while (b < numBins && value >= breaks[b]);

            observed[b - 1]++;
        }

        final double[] expected = new double[numBins];
        Arrays.fill(expected, (double) values.length / numBins);

        return new GTest().gTest(expected, observed);
    }
}