/*
 * Copyright (c) NM LTD.
 * https://nm.dev/
 * 
 * THIS SOFTWARE IS LICENSED, NOT SOLD.
 * 
 * YOU MAY USE THIS SOFTWARE ONLY AS DESCRIBED IN THE LICENSE.
 * IF YOU ARE NOT AWARE OF AND/OR DO NOT AGREE TO THE TERMS OF THE LICENSE,
 * DO NOT USE THIS SOFTWARE.
 * 
 * THE SOFTWARE IS PROVIDED "AS IS", WITH NO WARRANTY WHATSOEVER,
 * EITHER EXPRESS OR IMPLIED, INCLUDING, WITHOUT LIMITATION,
 * ANY WARRANTIES OF ACCURACY, ACCESSIBILITY, COMPLETENESS,
 * FITNESS FOR A PARTICULAR PURPOSE, MERCHANTABILITY, NON-INFRINGEMENT, 
 * TITLE AND USEFULNESS.
 * 
 * IN NO EVENT AND UNDER NO LEGAL THEORY,
 * WHETHER IN ACTION, CONTRACT, NEGLIGENCE, TORT, OR OTHERWISE,
 * SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR
 * ANY CLAIMS, DAMAGES OR OTHER LIABILITIES,
 * ARISING AS A RESULT OF USING OR OTHER DEALINGS IN THE SOFTWARE.
 */
package dev.nm.nmj;

import dev.nm.algebra.linear.matrix.doubles.Matrix;
import dev.nm.algebra.linear.matrix.doubles.matrixtype.dense.DenseMatrix;
import dev.nm.algebra.linear.vector.doubles.Vector;
import dev.nm.algebra.linear.vector.doubles.dense.DenseVector;
import dev.nm.analysis.function.rn2r1.univariate.AbstractUnivariateRealFunction;
import dev.nm.analysis.function.rn2r1.univariate.UnivariateRealFunction;
import dev.nm.analysis.function.special.gaussian.CumulativeNormalMarsaglia;
import dev.nm.analysis.function.special.gaussian.Gaussian;
import dev.nm.analysis.function.special.gaussian.StandardCumulativeNormal;
import dev.nm.interval.RealInterval;
import dev.nm.number.DoubleUtils;
import dev.nm.stat.descriptive.covariance.SampleCovariance;
import dev.nm.stat.descriptive.moment.Kurtosis;
import dev.nm.stat.descriptive.moment.Mean;
import dev.nm.stat.descriptive.moment.Skewness;
import dev.nm.stat.descriptive.moment.Variance;
import dev.nm.stat.distribution.univariate.BetaDistribution;
import dev.nm.stat.distribution.univariate.EmpiricalDistribution;
import dev.nm.stat.distribution.univariate.ExponentialDistribution;
import dev.nm.stat.distribution.univariate.GammaDistribution;
import dev.nm.stat.distribution.univariate.PoissonDistribution;
import dev.nm.stat.distribution.univariate.ProbabilityDistribution;
import dev.nm.stat.random.Estimator;
import dev.nm.stat.random.rng.multivariate.MultinomialRVG;
import dev.nm.stat.random.rng.multivariate.NormalRVG;
import dev.nm.stat.random.rng.multivariate.RandomVectorGenerator;
import dev.nm.stat.random.rng.multivariate.UniformDistributionOverBox;
import dev.nm.stat.random.rng.univariate.InverseTransformSampling;
import dev.nm.stat.random.rng.univariate.RandomLongGenerator;
import dev.nm.stat.random.rng.univariate.RandomNumberGenerator;
import dev.nm.stat.random.rng.univariate.beta.Cheng1978;
import dev.nm.stat.random.rng.univariate.beta.RandomBetaGenerator;
import dev.nm.stat.random.rng.univariate.exp.RandomExpGenerator;
import dev.nm.stat.random.rng.univariate.exp.Ziggurat2000Exp;
import dev.nm.stat.random.rng.univariate.gamma.KunduGupta2007;
import dev.nm.stat.random.rng.univariate.normal.NormalRNG;
import dev.nm.stat.random.rng.univariate.normal.RandomStandardNormalGenerator;
import dev.nm.stat.random.rng.univariate.normal.Zignor2005;
import dev.nm.stat.random.rng.univariate.poisson.Knuth1969;
import dev.nm.stat.random.rng.univariate.uniform.UniformRNG;
import dev.nm.stat.random.rng.univariate.uniform.linear.CompositeLinearCongruentialGenerator;
import dev.nm.stat.random.rng.univariate.uniform.linear.LEcuyer;
import dev.nm.stat.random.rng.univariate.uniform.linear.Lehmer;
import dev.nm.stat.random.rng.univariate.uniform.linear.LinearCongruentialGenerator;
import dev.nm.stat.random.rng.univariate.uniform.mersennetwister.MersenneTwister;
import dev.nm.stat.random.sampler.resampler.BootstrapEstimator;
import dev.nm.stat.random.sampler.resampler.bootstrap.CaseResamplingReplacement;
import dev.nm.stat.random.sampler.resampler.bootstrap.block.PattonPolitisWhite2009;
import dev.nm.stat.random.sampler.resampler.bootstrap.block.PattonPolitisWhite2009ForObject;
import dev.nm.stat.random.variancereduction.AntitheticVariates;
import dev.nm.stat.random.variancereduction.CommonRandomNumbers;
import dev.nm.stat.random.variancereduction.ControlVariates;
import dev.nm.stat.random.variancereduction.ImportanceSampling;
import dev.nm.stat.test.distribution.normality.ShapiroWilk;
import static java.lang.Math.PI;
import static java.lang.Math.sin;
import static java.lang.Math.sqrt;
import java.util.Arrays;

/**
 * Numerical Methods Using Java: For Data Science, Analysis, and Engineering
 *
 * @author haksunli
 * @see
 * https://www.amazon.com/Numerical-Methods-Using-Java-Engineering/dp/1484267966
 * https://nm.dev/
 */
public class Chapter13 {

    public static void main(String[] args) {
        System.out.println("Chapter 13 demos");

        Chapter13 chapter13 = new Chapter13();
        chapter13.lcgs();
        chapter13.MT19937();
        chapter13.normal_rng();
        chapter13.beta_rng();
        chapter13.gamma_rng();
        chapter13.Poisson_rng();
        chapter13.exponential_rng();
        chapter13.compute_Pi();
        chapter13.normal_rvg();
        chapter13.multinomial_rvg();
        chapter13.empirical_rng();
        chapter13.case_resampling_1();
        chapter13.case_resampling_2();
        chapter13.bootstrapping_methods();
        chapter13.crn();
        chapter13.antithetic_variates();
        chapter13.control_variates();
        chapter13.importance_sampling_1();
        chapter13.importance_sampling_2();
    }

    // from
    // https://www.scratchapixel.com/lessons/mathematics-physics-for-computer-graphics/monte-carlo-methods-in-practice/variance-reduction-methods
    private void importance_sampling_2() {
        RandomNumberGenerator rng = new UniformRNG();
        rng.seed(1234567892L);

        int N = 16;
        for (int n = 0; n < 10; ++n) {
            float sumUniform = 0, sumImportance = 0;
            for (int i = 0; i < N; ++i) {
                double r = rng.nextDouble();
                sumUniform += sin(r * PI * 0.5);
                double xi = sqrt(r) * PI * 0.5;
                sumImportance += sin(xi) / ((8 * xi) / (PI * PI));
            }
            sumUniform *= (PI * 0.5) / N;
            sumImportance *= 1.f / N;
            System.out.println(String.format("%f %f\n", sumUniform, sumImportance));
        }
    }

    private void importance_sampling_1() {
        UnivariateRealFunction h = new AbstractUnivariateRealFunction() {

            @Override
            public double evaluate(double x) {
                return x; // the identity function
            }
        };

        UnivariateRealFunction w = new AbstractUnivariateRealFunction() {
            private final Gaussian phi = new Gaussian();
            private final StandardCumulativeNormal N = new CumulativeNormalMarsaglia();
            private final double I = N.evaluate(1) - N.evaluate(0);

            @Override
            public double evaluate(double x) {
                double w = phi.evaluate(x) / I; // the weight
                return w;
            }
        };

        RandomNumberGenerator rng = new UniformRNG();
        rng.seed(1234567892L);

        ImportanceSampling is = new ImportanceSampling(h, w, rng);
        Estimator estimator = is.estimate(100000);
        System.out.println(
                String.format(
                        "mean = %f, variance = %f",
                        estimator.mean(),
                        estimator.variance()));
    }

    private void control_variates() {
        UnivariateRealFunction f
                = new AbstractUnivariateRealFunction() {

            @Override
            public double evaluate(double x) {
                double fx = 1. / (1. + x);
                return fx;
            }
        };

        UnivariateRealFunction g
                = new AbstractUnivariateRealFunction() {

            @Override
            public double evaluate(double x) {
                double gx = 1. + x;
                return gx;
            }
        };

        RandomLongGenerator uniform = new UniformRNG();
        uniform.seed(1234567891L);

        ControlVariates cv
                = new ControlVariates(f, g, 1.5, -0.4773, uniform);
        ControlVariates.Estimator estimator = cv.estimate(1500);
        System.out.println(
                String.format(
                        "mean = %f, variance = %f, b = %f",
                        estimator.mean(),
                        estimator.variance(),
                        estimator.b()));
    }

    private void antithetic_variates() {
        UnivariateRealFunction f
                = new AbstractUnivariateRealFunction() {

            @Override
            public double evaluate(double x) {
                double fx = 1. / (1. + x);
                return fx;
            }
        };

        RandomLongGenerator uniform = new UniformRNG();
        uniform.seed(1234567894L);

        AntitheticVariates av
                = new AntitheticVariates(
                        f,
                        uniform,
                        AntitheticVariates.REFLECTION);
        Estimator estimator = av.estimate(1500);
        System.out.println(
                String.format(
                        "mean = %f, variance = %f",
                        estimator.mean(),
                        estimator.variance()));
    }

    private void crn() {
        final UnivariateRealFunction f
                = new AbstractUnivariateRealFunction() {

            @Override
            public double evaluate(double x) {
                double fx = 2. - Math.sin(x) / x;
                return fx;
            }
        };

        final UnivariateRealFunction g
                = new AbstractUnivariateRealFunction() {

            @Override
            public double evaluate(double x) {
                double gx = Math.exp(x * x) - 0.5;
                return gx;
            }
        };

        RandomLongGenerator X1 = new UniformRNG();
        X1.seed(1234567890L);

        CommonRandomNumbers crn0
                = new CommonRandomNumbers(
                        f,
                        g,
                        X1,
                        new AbstractUnivariateRealFunction() { // another independent uniform RNG
                    final RandomLongGenerator X2 = new UniformRNG();

                    {
                        X2.seed(246890123L);
                    }

                    @Override
                    public double evaluate(double x) {
                        return X2.nextDouble();
                    }
                });
        Estimator estimator0 = crn0.estimate(100_000);
        System.out.println(
                String.format("d = %f, variance = %f",
                        estimator0.mean(),
                        estimator0.variance()));

        CommonRandomNumbers crn1
                = new CommonRandomNumbers(f, g, X1); // use X1 for both f and g
        Estimator estimator1 = crn1.estimate(100_000);
        System.out.println(
                String.format("d = %f, variance = %f",
                        estimator1.mean(),
                        estimator1.variance()));
    }

    /**
     * Constructs a dependent sequence (consisting of 0 and 1) by retaining the
     * last value with probability <i>q</i> while changing the last value with
     * probability <i>1-q</i>.
     * <p/>
     * The simple bootstrapping method {@linkplain CaseResamplingReplacement}
     * will severely overestimate the occurrences of certain pattern, while
     * block bootstrapping method {@linkplain BlockBootstrap} gives a good
     * estimation of the occurrences in the original sample. All estimators over
     * estimate.
     */
    private void bootstrapping_methods() {
        final int N = 10000;
        final double q = 0.70; // the probability of retaining last value

        UniformRNG uniformRNG = new UniformRNG();
        uniformRNG.seed(1234567890L);

        // generate a randome series of 0s and 1s with serial correlation
        final double[] sample = new double[N];
        sample[0] = uniformRNG.nextDouble() > 0.5 ? 1 : 0;
        for (int i = 1; i < N; ++i) {
            sample[i] = uniformRNG.nextDouble() < q ? sample[i - 1] : 1 - sample[i - 1];
        }

        // simple case resampling with replacement method
        CaseResamplingReplacement simpleBoot
                = new CaseResamplingReplacement(sample, uniformRNG);
        Mean countInSimpleBootstrap = new Mean();

        RandomNumberGenerator rlg = new Ziggurat2000Exp();
        rlg.seed(1234567890L);

        // Patton-Politis-White method using stationary blocks
        PattonPolitisWhite2009 stationaryBlock
                = new PattonPolitisWhite2009(
                        sample,
                        PattonPolitisWhite2009ForObject.Type.STATIONARY,
                        uniformRNG,
                        rlg);
        Mean countInStationaryBlockBootstrap = new Mean();

        // Patton-Politis-White method using circular blocks
        PattonPolitisWhite2009 circularBlock
                = new PattonPolitisWhite2009(
                        sample,
                        PattonPolitisWhite2009ForObject.Type.CIRCULAR,
                        uniformRNG,
                        rlg);
        Mean countInCircularBlockBootstrap = new Mean();

        // change this line to use a different pattern
        final double[] pattern = new double[]{1, 0, 1, 0, 1};

        final int B = 10000;
        for (int i = 0; i < B; ++i) {
            // count the number of occurrences for the pattern in the series
            int numberOfMatches = match(simpleBoot.newResample(), pattern);
            countInSimpleBootstrap.addData(numberOfMatches);

            // count the number of occurrences for the pattern in the series
            numberOfMatches = match(stationaryBlock.newResample(), pattern);
            countInStationaryBlockBootstrap.addData(numberOfMatches);

            // count the number of occurrences for the pattern in the series
            numberOfMatches = match(circularBlock.newResample(), pattern);
            countInCircularBlockBootstrap.addData(numberOfMatches);
        }

        // compare the numbers of occurrences of the pattern using different bootstrap methods
        int countInSample = match(sample, pattern);
        System.out.println("matched patterns in sample: " + countInSample);
        System.out.println("matched patterns in simple bootstrap: " + countInSimpleBootstrap.value());
        System.out.println("matched patterns in stationary block bootstrap: " + countInStationaryBlockBootstrap.value());
        System.out.println("matched patterns in circular block bootstrap: " + countInCircularBlockBootstrap.value());
    }

    private static int match(double[] seq, double[] pattern) {
        int count = 0;
        for (int i = 0; i < seq.length - pattern.length; ++i) {
            if (seq[i] == pattern[0]) {
                double[] trunc = Arrays.copyOfRange(seq, i, i + pattern.length);
                if (DoubleUtils.equal(trunc, pattern, 1e-7)) {
                    count++;
                }
            }
        }
        return count;
    }

    private void case_resampling_1() {
        // sample from true population
        double[] sample = new double[]{150., 155., 160., 165., 170.};
        CaseResamplingReplacement boot = new CaseResamplingReplacement(sample);
        boot.seed(1234567890L);

        int B = 1000;
        double[] means = new double[B];
        for (int i = 0; i < B; ++i) {
            double[] resample = boot.newResample();
            means[i] = new Mean(resample).value();
        }

        // estimator of population mean
        double mean = new Mean(means).value();
        // variance of estimator; limited by sample size (regardless of how big B is)
        double var = new Variance(means).value();

        System.out.println(
                String.format("mean = %f, variance of the estimated mean = %f",
                        mean,
                        var));
    }

    private void case_resampling_2() {
        // sample from true population
        double[] sample = new double[]{150., 155., 160., 165., 170.};
        CaseResamplingReplacement boot = new CaseResamplingReplacement(sample);
        boot.seed(1234567890L);

        int B = 1000;
        BootstrapEstimator estimator
                = new BootstrapEstimator(boot, () -> new Mean(), B);

        System.out.println(
                String.format("mean = %f, variance of the estimated mean = %f",
                        estimator.value(),
                        estimator.variance()));
    }

    private void empirical_rng() {
        // we first generate some samples from standard normal distribution
        RandomLongGenerator uniform = new MersenneTwister();
        uniform.seed(1234567890L);

        RandomStandardNormalGenerator rng1 = new Zignor2005(uniform); // mean = 0, stdev = 1
        int N = 1000;
        double[] x1 = new double[N];
        for (int i = 0; i < N; ++i) {
            x1[i] = rng1.nextDouble();
        }

        // compute the empirical distribution function from the sample data
        EmpiricalDistribution dist2 = new EmpiricalDistribution(x1);

        // construct an RNG using inverse transform sampling method
        InverseTransformSampling rng2 = new InverseTransformSampling(dist2);

        // generate some random variates from the RNG
        double[] x2 = new double[N];
        for (int i = 0; i < N; ++i) {
            x2[i] = rng2.nextDouble();
        }

        // check the properties of the random variates
        Variance var = new Variance(x2);
        double mean = var.mean();
        double stdev = var.standardDeviation();
        System.out.println(String.format("mean = %f, standard deviation = %f", mean, stdev));

        // check if the samples are normally distributed
        ShapiroWilk test = new ShapiroWilk(x2);
        System.out.println(String.format("ShapiroWilk statistics = %f, pValue = %f", test.statistics(), test.pValue()));
    }

    private void multinomial_rvg() {
        MultinomialRVG rvg
                = new MultinomialRVG(100_000, new double[]{0.7, 0.3}); // bin0 is 70% chance, bin1 30% chance
        double[] bin = rvg.nextVector();

        double total = 0;
        for (int i = 0; i < bin.length; ++i) {
            total += bin[i];
        }

        double bin0 = bin[0] / total; // bin0 percentage
        double bin1 = bin[1] / total; // bin0 percentage

        System.out.println(String.format("bin0 %% = %f, bin1 %% = %f", bin0, bin1));
    }

    private void normal_rvg() {
        // mean
        Vector mu = new DenseVector(new double[]{-2., 2.});
        // covariance matrix
        Matrix sigma = new DenseMatrix(new double[][]{
            {1., 0.5},
            {0.5, 1.}
        });
        NormalRVG rvg = new NormalRVG(mu, sigma);
        rvg.seed(1234567890L);

        final int size = 10_000;
        double[][] x = new double[size][];

        Mean mean1 = new Mean();
        Mean mean2 = new Mean();
        for (int i = 0; i < size; ++i) {
            double[] v = rvg.nextVector();
            mean1.addData(v[0]);
            mean2.addData(v[1]);

            x[i] = v;
        }

        System.out.println(String.format("mean of X_1 = %f", mean1.value()));
        System.out.println(String.format("mean of X_2 = %f", mean2.value()));

        Matrix X = new DenseMatrix(x);
        SampleCovariance cov = new SampleCovariance(X);
        System.out.println(String.format("sample covariance = %s", cov.toString()));
    }

    private void compute_Pi() {
        final int N = 1_000_000;

        RandomVectorGenerator rvg
                = new UniformDistributionOverBox(
                        new RealInterval(-1., 1.), // a unit square box
                        new RealInterval(-1., 1.));

        int N0 = 0;
        for (int i = 0; i < N; i++) {
            double[] xy = rvg.nextVector();
            double x = xy[0], y = xy[1];
            if (x * x + y * y <= 1.) { // check if the dot is inside a circle
                N0++;
            }
        }
        double pi = 4. * N0 / N;
        System.out.println("pi = " + pi);
    }

    private void exponential_rng() {
        int size = 500_000;

        RandomExpGenerator rng = new Ziggurat2000Exp();
        rng.seed(634641070L);

        double[] x = new double[size];
        for (int i = 0; i < size; ++i) {
            x[i] = rng.nextDouble();
        }

        // compute the sample statistics
        Mean mean = new Mean(x);
        Variance var = new Variance(x);
        Skewness skew = new Skewness(x);
        Kurtosis kurtosis = new Kurtosis(x);

        // compute the theoretial statistics
        ProbabilityDistribution dist = new ExponentialDistribution();

        // compute the theoretial statistics
        printStats(dist, mean, var, skew, kurtosis);
    }

    private void Poisson_rng() {
        final int N = 10_000;
        double lambda = 1;

        RandomNumberGenerator rng = new Knuth1969(lambda);
        rng.seed(123456789L);

        double[] x = new double[N];
        for (int i = 0; i < N; ++i) {
            x[i] = rng.nextDouble();
        }

        // compute the sample statistics
        Mean mean = new Mean(x);
        Variance var = new Variance(x);
        Skewness skew = new Skewness(x);
        Kurtosis kurtosis = new Kurtosis(x);

        // compute the theoretial statistics
        PoissonDistribution dist = new PoissonDistribution(lambda);

        // compute the theoretial statistics
        printStats(dist, mean, var, skew, kurtosis);
    }

    private void gamma_rng() {
        final int size = 1_000_000;

        final double k = 0.1;
        final double theta = 1;

        KunduGupta2007 rng = new KunduGupta2007(k, theta, new UniformRNG());
        rng.seed(1234567895L);

        double[] x = new double[size];
        for (int i = 0; i < size; ++i) {
            x[i] = rng.nextDouble();
        }

        // compute the sample statistics
        Mean mean = new Mean(x);
        Variance var = new Variance(x);
        Skewness skew = new Skewness(x);
        Kurtosis kurtosis = new Kurtosis(x);

        // compute the theoretial statistics
        ProbabilityDistribution dist = new GammaDistribution(k, theta);

        // compute the theoretial statistics
        printStats(dist, mean, var, skew, kurtosis);
    }

    private void beta_rng() {
        final int size = 1_000_000;

        final double alpha = 0.1;
        final double beta = 0.2;

        RandomBetaGenerator rng = new Cheng1978(alpha, beta, new UniformRNG());
        rng.seed(1234567890L);

        double[] x = new double[size];
        for (int i = 0; i < size; ++i) {
            x[i] = rng.nextDouble();
        }

        // compute the sample statistics
        Mean mean = new Mean(x);
        Variance var = new Variance(x);
        Skewness skew = new Skewness(x);
        Kurtosis kurtosis = new Kurtosis(x);

        // compute the theoretial statistics
        ProbabilityDistribution dist = new BetaDistribution(alpha, beta);

        // compare sample vs theoretical statistics
        printStats(dist, mean, var, skew, kurtosis);
    }

    private void normal_rng() {
        RandomLongGenerator uniform = new MersenneTwister();
        uniform.seed(1234567890L);

        RandomStandardNormalGenerator rng1 = new Zignor2005(uniform); // mean = 0, stdev = 1
        int N = 1000;
        double[] arr1 = new double[N];
        for (int i = 0; i < N; ++i) {
            arr1[i] = rng1.nextDouble();
        }

        // check the statistics of the random samples
        Variance var1 = new Variance(arr1);
        System.out.println(
                String.format(
                        "mean = %f, stdev = %f",
                        var1.mean(),
                        var1.standardDeviation()));

        NormalRNG rng2 = new NormalRNG(1., 2., rng1); // mean = 1, stdev = 2
        double[] arr2 = new double[N];
        for (int i = 0; i < N; ++i) {
            arr2[i] = rng2.nextDouble();
        }

        // check the statistics of the random samples
        Variance var2 = new Variance(arr2);
        System.out.println(
                String.format(
                        "mean = %f, stdev = %f",
                        var2.mean(),
                        var2.standardDeviation()));
    }

    private void MT19937() {
        RandomLongGenerator rng = new MersenneTwister();

        long startTime = System.nanoTime();
        int N = 1_000_000;
        for (int i = 0; i < N; ++i) {
            rng.nextDouble();
        }
        long endTime = System.nanoTime();

        long duration = (endTime - startTime);
        double ms = (double) duration / 1_000_000.; // divide by 1000000 to get milliseconds
        System.out.println(String.format("took MT19937 %f milliseconds to generate %d random numbers", ms, N));
    }

    private void lcgs() {
        System.out.println("generate randome numbers using an Lehmer RNG:");
        RandomLongGenerator rng1 = new Lehmer();
        rng1.seed(1234567890L);
        generateIntAndPrint(rng1, 10);
        double[] arr = generate(rng1, 10);
        print(arr);

        System.out.println("generate randome numbers using an LEcuyer RNG:");
        RandomLongGenerator rng2 = new LEcuyer();
        rng2.seed(1234567890L);
        generateIntAndPrint(rng2, 10);
        arr = generate(rng2, 10);
        print(arr);

        System.out.println("generate randome numbers using a composite LCG:");
        RandomLongGenerator rng3
                = new CompositeLinearCongruentialGenerator(
                        new LinearCongruentialGenerator[]{
                            (LinearCongruentialGenerator) rng1,
                            (LinearCongruentialGenerator) rng2
                        }
                );
        rng3.seed(1234567890L);
        generateIntAndPrint(rng3, 10);
        arr = generate(rng3, 10);
        print(arr);
    }

    private static double[] generate(RandomNumberGenerator rng, int n) {
        double[] arr = new double[n];
        for (int i = 0; i < n; i++) {
            arr[i] = rng.nextDouble();
        }
        return arr;
    }

    private static void print(double[] arr) {
        System.out.println(Arrays.toString(arr));
    }

    private static void generateIntAndPrint(RandomNumberGenerator rng, int n) {
        double[] randomNumbers = new double[n];
        for (int i = 0; i < n; i++) {
            randomNumbers[i] = rng.nextDouble();
        }
        System.out.println(Arrays.toString(randomNumbers));
    }

    private void printStats(
            ProbabilityDistribution dist,
            Mean mean,
            Variance var,
            Skewness skew,
            Kurtosis kurtosis
    ) {
        System.out.println(
                String.format("theoretical mean = %f, sample mean = %f",
                        dist.mean(),
                        mean.value()));
        System.out.println(
                String.format("theoretical var = %f, sample var = %f",
                        dist.variance(),
                        var.value()));
        System.out.println(
                String.format("theoretical skew = %f, sample skew = %f",
                        dist.skew(),
                        skew.value()));
        System.out.println(
                String.format("theoretical kurtosis = %f, sample kurtosis = %f",
                        dist.kurtosis(),
                        kurtosis.value()));
    }
}