Updated to use factory constructors for Statistics distributions
diff --git a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/distribution/EmpiricalDistribution.java b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/distribution/EmpiricalDistribution.java
index cca2c19..8058ea5 100644
--- a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/distribution/EmpiricalDistribution.java
+++ b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/distribution/EmpiricalDistribution.java
@@ -551,8 +551,8 @@
                 stats.getVariance() == 0) {
                 return new ConstantContinuousDistribution(stats.getMean());
             } else {
-                return new NormalDistribution(stats.getMean(),
-                                              stats.getStandardDeviation());
+                return NormalDistribution.of(stats.getMean(),
+                                             stats.getStandardDeviation());
             }
         };
     }
diff --git a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/distribution/MultivariateNormalDistribution.java b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/distribution/MultivariateNormalDistribution.java
index 6a72dbb..b216816 100644
--- a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/distribution/MultivariateNormalDistribution.java
+++ b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/distribution/MultivariateNormalDistribution.java
@@ -181,7 +181,7 @@
     public MultivariateRealDistribution.Sampler createSampler(final UniformRandomProvider rng) {
         return new MultivariateRealDistribution.Sampler() {
             /** Normal distribution. */
-            private final ContinuousDistribution.Sampler gauss = new NormalDistribution(0, 1).createSampler(rng);
+            private final ContinuousDistribution.Sampler gauss = NormalDistribution.of(0, 1).createSampler(rng);
 
             /** {@inheritDoc} */
             @Override
diff --git a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/CMAESOptimizer.java b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/CMAESOptimizer.java
index 83b036e..eb92f7b 100644
--- a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/CMAESOptimizer.java
+++ b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/optim/nonlinear/scalar/noderiv/CMAESOptimizer.java
@@ -246,7 +246,7 @@
         this.isActiveCMA = isActiveCMA;
         this.diagonalOnly = diagonalOnly;
         this.checkFeasableCount = Math.max(0, checkFeasableCount);
-        this.random = new NormalDistribution(0, 1).createSampler(rng);
+        this.random = NormalDistribution.of(0, 1).createSampler(rng);
         this.generateStatistics = generateStatistics;
     }
 
diff --git a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/correlation/PearsonsCorrelation.java b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/correlation/PearsonsCorrelation.java
index 34c2f6d..eda73a6 100644
--- a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/correlation/PearsonsCorrelation.java
+++ b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/correlation/PearsonsCorrelation.java
@@ -192,7 +192,7 @@
      * @throws NullPointerException if this instance was created with no data.
      */
     public RealMatrix getCorrelationPValues() {
-        TDistribution tDistribution = new TDistribution(nObs - 2);
+        TDistribution tDistribution = TDistribution.of(nObs - 2);
         int nVars = correlationMatrix.getColumnDimension();
         double[][] out = new double[nVars][nVars];
         for (int i = 0; i < nVars; i++) {
diff --git a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/inference/BinomialTest.java b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/inference/BinomialTest.java
index 011b7aa..80650c9 100644
--- a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/inference/BinomialTest.java
+++ b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/inference/BinomialTest.java
@@ -119,7 +119,7 @@
             throw new NullArgumentException();
         }
 
-        final BinomialDistribution distribution = new BinomialDistribution(numberOfTrials, probability);
+        final BinomialDistribution distribution = BinomialDistribution.of(numberOfTrials, probability);
         switch (alternativeHypothesis) {
         case GREATER_THAN:
             return 1 - distribution.cumulativeProbability(numberOfSuccesses - 1);
diff --git a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/inference/ChiSquareTest.java b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/inference/ChiSquareTest.java
index 92f43a0..e695d3f 100644
--- a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/inference/ChiSquareTest.java
+++ b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/inference/ChiSquareTest.java
@@ -157,7 +157,7 @@
 
         // pass a null rng to avoid unneeded overhead as we will not sample from this distribution
         final ChiSquaredDistribution distribution =
-            new ChiSquaredDistribution(expected.length - 1.0);
+            ChiSquaredDistribution.of(expected.length - 1.0);
         return 1.0 - distribution.cumulativeProbability(chiSquare(expected, observed));
     }
 
@@ -332,7 +332,7 @@
         checkArray(counts);
         double df = ((double) counts.length -1) * ((double) counts[0].length - 1);
         // pass a null rng to avoid unneeded overhead as we will not sample from this distribution
-        final ChiSquaredDistribution distribution = new ChiSquaredDistribution(df);
+        final ChiSquaredDistribution distribution = ChiSquaredDistribution.of(df);
         return 1 - distribution.cumulativeProbability(chiSquare(counts));
 
     }
@@ -536,7 +536,7 @@
 
         // pass a null rng to avoid unneeded overhead as we will not sample from this distribution
         final ChiSquaredDistribution distribution =
-                new ChiSquaredDistribution((double) observed1.length - 1);
+                ChiSquaredDistribution.of((double) observed1.length - 1);
         return 1 - distribution.cumulativeProbability(
                 chiSquareDataSetsComparison(observed1, observed2));
 
diff --git a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/inference/GTest.java b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/inference/GTest.java
index baaefb3..6dddb4a 100644
--- a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/inference/GTest.java
+++ b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/inference/GTest.java
@@ -154,7 +154,7 @@
 
         // pass a null rng to avoid unneeded overhead as we will not sample from this distribution
         final ChiSquaredDistribution distribution =
-                new ChiSquaredDistribution(expected.length - 1.0);
+                ChiSquaredDistribution.of(expected.length - 1.0);
         return 1.0 - distribution.cumulativeProbability(g(expected, observed));
     }
 
@@ -185,7 +185,7 @@
 
         // pass a null rng to avoid unneeded overhead as we will not sample from this distribution
         final ChiSquaredDistribution distribution =
-                new ChiSquaredDistribution(expected.length - 2.0);
+                ChiSquaredDistribution.of(expected.length - 2.0);
         return 1.0 - distribution.cumulativeProbability(g(expected, observed));
     }
 
@@ -475,7 +475,7 @@
 
         // pass a null rng to avoid unneeded overhead as we will not sample from this distribution
         final ChiSquaredDistribution distribution =
-                new ChiSquaredDistribution((double) observed1.length - 1);
+                ChiSquaredDistribution.of((double) observed1.length - 1);
         return 1 - distribution.cumulativeProbability(
                 gDataSetsComparison(observed1, observed2));
     }
diff --git a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/inference/MannWhitneyUTest.java b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/inference/MannWhitneyUTest.java
index be14563..2f1cdd3 100644
--- a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/inference/MannWhitneyUTest.java
+++ b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/inference/MannWhitneyUTest.java
@@ -181,7 +181,7 @@
 
         // No try-catch or advertised exception because args are valid
         // pass a null rng to avoid unneeded overhead as we will not sample from this distribution
-        final NormalDistribution standardNormal = new NormalDistribution(0, 1);
+        final NormalDistribution standardNormal = NormalDistribution.of(0, 1);
 
         return 2 * standardNormal.cumulativeProbability(z);
     }
diff --git a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/inference/OneWayAnova.java b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/inference/OneWayAnova.java
index e76cb20..d477427 100644
--- a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/inference/OneWayAnova.java
+++ b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/inference/OneWayAnova.java
@@ -126,7 +126,7 @@
         final AnovaStats a = anovaStats(categoryData);
         // No try-catch or advertised exception because args are valid
         // pass a null rng to avoid unneeded overhead as we will not sample from this distribution
-        final FDistribution fdist = new FDistribution(a.dfbg, a.dfwg);
+        final FDistribution fdist = FDistribution.of(a.dfbg, a.dfwg);
         return 1.0 - fdist.cumulativeProbability(a.f);
 
     }
@@ -168,7 +168,7 @@
 
         final AnovaStats a = anovaStats(categoryData, allowOneElementData);
         // pass a null rng to avoid unneeded overhead as we will not sample from this distribution
-        final FDistribution fdist = new FDistribution(a.dfbg, a.dfwg);
+        final FDistribution fdist = FDistribution.of(a.dfbg, a.dfwg);
         return 1.0 - fdist.cumulativeProbability(a.f);
 
     }
diff --git a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/inference/TTest.java b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/inference/TTest.java
index 9854538..be4dc90 100644
--- a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/inference/TTest.java
+++ b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/inference/TTest.java
@@ -1057,7 +1057,7 @@
 
         final double t = AccurateMath.abs(t(m, mu, v, n));
         // pass a null rng to avoid unneeded overhead as we will not sample from this distribution
-        final TDistribution distribution = new TDistribution(n - 1);
+        final TDistribution distribution = TDistribution.of(n - 1);
         return 2.0 * distribution.cumulativeProbability(-t);
 
     }
@@ -1087,7 +1087,7 @@
         final double t = AccurateMath.abs(t(m1, m2, v1, v2, n1, n2));
         final double degreesOfFreedom = df(v1, v2, n1, n2);
         // pass a null rng to avoid unneeded overhead as we will not sample from this distribution
-        final TDistribution distribution = new TDistribution(degreesOfFreedom);
+        final TDistribution distribution = TDistribution.of(degreesOfFreedom);
         return 2.0 * distribution.cumulativeProbability(-t);
 
     }
@@ -1117,7 +1117,7 @@
         final double t = AccurateMath.abs(homoscedasticT(m1, m2, v1, v2, n1, n2));
         final double degreesOfFreedom = n1 + n2 - 2;
         // pass a null rng to avoid unneeded overhead as we will not sample from this distribution
-        final TDistribution distribution = new TDistribution(degreesOfFreedom);
+        final TDistribution distribution = TDistribution.of(degreesOfFreedom);
         return 2.0 * distribution.cumulativeProbability(-t);
 
     }
diff --git a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/inference/WilcoxonSignedRankTest.java b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/inference/WilcoxonSignedRankTest.java
index 8d7141d..0d1fdd4 100644
--- a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/inference/WilcoxonSignedRankTest.java
+++ b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/inference/WilcoxonSignedRankTest.java
@@ -253,7 +253,7 @@
 
         // No try-catch or advertised exception because args are valid
         // pass a null rng to avoid unneeded overhead as we will not sample from this distribution
-        final NormalDistribution standardNormal = new NormalDistribution(0, 1);
+        final NormalDistribution standardNormal = NormalDistribution.of(0, 1);
 
         return 2*standardNormal.cumulativeProbability(z);
     }
diff --git a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/interval/AgrestiCoullInterval.java b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/interval/AgrestiCoullInterval.java
index 3ddb4cb..2e60415 100644
--- a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/interval/AgrestiCoullInterval.java
+++ b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/interval/AgrestiCoullInterval.java
@@ -34,7 +34,7 @@
     public ConfidenceInterval createInterval(int numberOfTrials, int numberOfSuccesses, double confidenceLevel) {
         IntervalUtils.checkParameters(numberOfTrials, numberOfSuccesses, confidenceLevel);
         final double alpha = (1.0 - confidenceLevel) / 2;
-        final NormalDistribution normalDistribution = new NormalDistribution(0, 1);
+        final NormalDistribution normalDistribution = NormalDistribution.of(0, 1);
         final double z = normalDistribution.inverseCumulativeProbability(1 - alpha);
         final double zSquared = AccurateMath.pow(z, 2);
         final double modifiedNumberOfTrials = numberOfTrials + zSquared;
diff --git a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/interval/ClopperPearsonInterval.java b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/interval/ClopperPearsonInterval.java
index 20f4913..a28522c 100644
--- a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/interval/ClopperPearsonInterval.java
+++ b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/interval/ClopperPearsonInterval.java
@@ -40,16 +40,16 @@
         final double alpha = 0.5 * (1 - confidenceLevel);
 
         if (numberOfSuccesses > 0) {
-            final FDistribution distributionLowerBound = new FDistribution(2.0 * (numberOfTrials - numberOfSuccesses + 1),
-                                                                           2.0 * numberOfSuccesses);
+            final FDistribution distributionLowerBound = FDistribution.of(2.0 * (numberOfTrials - numberOfSuccesses + 1),
+                                                                          2.0 * numberOfSuccesses);
             final double fValueLowerBound = distributionLowerBound.inverseCumulativeProbability(1 - alpha);
             lowerBound = numberOfSuccesses /
                 (numberOfSuccesses + (numberOfTrials - numberOfSuccesses + 1) * fValueLowerBound);
         }
 
         if (numberOfSuccesses < numberOfTrials) {
-            final FDistribution distributionUpperBound = new FDistribution(2.0 * (numberOfSuccesses + 1),
-                                                                           2.0 * (numberOfTrials - numberOfSuccesses));
+            final FDistribution distributionUpperBound = FDistribution.of(2.0 * (numberOfSuccesses + 1),
+                                                                          2.0 * (numberOfTrials - numberOfSuccesses));
             final double fValueUpperBound = distributionUpperBound.inverseCumulativeProbability(1 - alpha);
             upperBound = (numberOfSuccesses + 1) * fValueUpperBound /
                 (numberOfTrials - numberOfSuccesses + (numberOfSuccesses + 1) * fValueUpperBound);
diff --git a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/interval/NormalApproximationInterval.java b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/interval/NormalApproximationInterval.java
index 75533e4..c9e4832 100644
--- a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/interval/NormalApproximationInterval.java
+++ b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/interval/NormalApproximationInterval.java
@@ -36,7 +36,7 @@
         IntervalUtils.checkParameters(numberOfTrials, numberOfSuccesses, confidenceLevel);
         final double mean = (double) numberOfSuccesses / (double) numberOfTrials;
         final double alpha = (1.0 - confidenceLevel) / 2;
-        final NormalDistribution normalDistribution = new NormalDistribution(0, 1);
+        final NormalDistribution normalDistribution = NormalDistribution.of(0, 1);
         final double difference = normalDistribution.inverseCumulativeProbability(1 - alpha) *
                                   AccurateMath.sqrt(1.0 / numberOfTrials * mean * (1 - mean));
         return new ConfidenceInterval(mean - difference, mean + difference, confidenceLevel);
diff --git a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/interval/WilsonScoreInterval.java b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/interval/WilsonScoreInterval.java
index 713f46e..0fa89ec 100644
--- a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/interval/WilsonScoreInterval.java
+++ b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/interval/WilsonScoreInterval.java
@@ -32,7 +32,7 @@
     public ConfidenceInterval createInterval(int numberOfTrials, int numberOfSuccesses, double confidenceLevel) {
         IntervalUtils.checkParameters(numberOfTrials, numberOfSuccesses, confidenceLevel);
         final double alpha = (1 - confidenceLevel) / 2;
-        final NormalDistribution normalDistribution = new NormalDistribution(0, 1);
+        final NormalDistribution normalDistribution = NormalDistribution.of(0, 1);
         final double z = normalDistribution.inverseCumulativeProbability(1 - alpha);
         final double zSquared = z * z;
         final double oneOverNumTrials = 1d / numberOfTrials;
diff --git a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/regression/SimpleRegression.java b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/regression/SimpleRegression.java
index c4b18dd..981f14d 100644
--- a/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/regression/SimpleRegression.java
+++ b/commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/stat/regression/SimpleRegression.java
@@ -698,7 +698,7 @@
                                           alpha, 0, 1);
         }
         // No advertised NotStrictlyPositiveException here - will return NaN above
-        TDistribution distribution = new TDistribution(n - 2d);
+        TDistribution distribution = TDistribution.of(n - 2d);
         return getSlopeStdErr() *
             distribution.inverseCumulativeProbability(1d - alpha / 2d);
     }
@@ -730,7 +730,7 @@
             return Double.NaN;
         }
         // No advertised NotStrictlyPositiveException here - will return NaN above
-        TDistribution distribution = new TDistribution(n - 2d);
+        TDistribution distribution = TDistribution.of(n - 2d);
         return 2d * (1.0 - distribution.cumulativeProbability(
                     AccurateMath.abs(getSlope()) / getSlopeStdErr()));
     }