Minor bug fixes in hypothesis tests, MGBT, binomial. Improving convergence in mixture and GMM expectation-maximization method. Made degrees of freedom in student-t and noncentral-t doubles.

HadenSmith · HadenSmith · commit ae06c340ac8c · 2026-03-28T06:54:11.000-06:00
diff --git a/Numerics/Data/Statistics/HypothesisTests.cs b/Numerics/Data/Statistics/HypothesisTests.cs
@@ -28,13 +28,14 @@
 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */
 
-using System;
-using System.Collections.Generic;
-using System.Linq;
 using Numerics.Distributions;
 using Numerics.MachineLearning;
 using Numerics.Mathematics.LinearAlgebra;
 using Numerics.Mathematics.SpecialFunctions;
+using System;
+using System.Collections.Generic;
+using System.Diagnostics;
+using System.Linq;
 
 namespace Numerics.Data.Statistics
 {
@@ -111,7 +112,7 @@ public static double UnequalVarianceTtest(IList<double> sample1, IList<double> s
             int n2 = sample2.Count;
             double t = Math.Abs((ave1 - ave2) / Math.Sqrt(var1 / n1 + var2 / n2));
             double df = Tools.Sqr(var1 / n1 + var2 / n2) / (Tools.Sqr(var1 / n1) / (n1 - 1) + Tools.Sqr(var2 / n2) / (n2 - 1));
-            var tdist = new StudentT((int)df);
+            var tdist = new StudentT(df);
             return (1d - tdist.CDF(t)) * 2d;
         }
 
@@ -134,7 +135,7 @@ public static double PairedTtest(IList<double> sample1, IList<double> sample2)
             cov /= (df = n - 1);
             sd = Math.Sqrt((var1 + var2 - 2.0 * cov) / n);
             double t = Math.Abs((ave1 - ave2) / sd);
-            var tdist = new StudentT((int)df);
+            var tdist = new StudentT(df);
             return (1d - tdist.CDF(t)) * 2d;
         }
 
@@ -348,29 +349,36 @@ public static double UnimodalityTest(IList<double> sample)
             int n = sample.Count;
             if (n < 10) throw new ArgumentException("The sample size must be greater than or equal to 10.");
 
-            // Fit 1-component GMM (unimodal)
-            var gmm1 = new GaussianMixtureModel(sample.ToArray(), 1);
-            gmm1.Train(12345, true);
-            var logLH1 = gmm1.LogLikelihood;
+            try
+            {
+                // Fit 1-component GMM (unimodal)
+                var gmm1 = new GaussianMixtureModel(sample.ToArray(), 1);
+                gmm1.Train(12345, true);
+                var logLH1 = gmm1.LogLikelihood;
 
-            // Fit 2-component GMM (potentially bimodal)
-            var gmm2 = new GaussianMixtureModel(sample.ToArray(), 2);
-            gmm2.Train(12345, true);
-            var logLH2 = gmm2.LogLikelihood;
+                // Fit 2-component GMM (potentially bimodal)
+                var gmm2 = new GaussianMixtureModel(sample.ToArray(), 2);
+                gmm2.Train(12345, true);
+                var logLH2 = gmm2.LogLikelihood;
 
-            // Compute test statistic: Likelihood Ratio Test
-            double testStat = 2 * (logLH2 - logLH1);
-            int df = 3;
+                // Compute test statistic: Likelihood Ratio Test
+                double testStat = 2 * (logLH2 - logLH1);
+                int df = 3;
 
-            // Compute p-value from chi-square distribution
-            var chiSquared = new ChiSquared(df);
-            double pval = 1.0 - chiSquared.CDF(testStat);
+                // Compute p-value from chi-square distribution
+                var chiSquared = new ChiSquared(df);
+                double pval = 1.0 - chiSquared.CDF(testStat);
 
-            return pval;
+                return pval;
+            }
+            catch (Exception ex) 
+            {
+                Debug.WriteLine($"Error in hypothesis testing: {ex.Message}");
+                return double.NaN;
+            }
 
         }
 
-
     }
 
 
diff --git a/Numerics/Data/Statistics/MultipleGrubbsBeckTest.cs b/Numerics/Data/Statistics/MultipleGrubbsBeckTest.cs
@@ -233,8 +233,19 @@ private static double FGGB(double PZR, int N, int R, double ETA)
             df = 2.0d * alpha;
             ncp = (MuMP - ZR) / Math.Sqrt(VarMP);
             q = -Math.Sqrt(MuS2 / VarMP) * EtaP;
-            var NCTDist = new NoncentralT(df, ncp);
-            ANS = 1.0d - NCTDist.CDF(q);
+            // Match Fortran FP_TNC_CDF: use normal approximation for df > 20
+            double cdfResult;
+            if (df > 20.0d)
+            {
+                double Z = (q * (1.0d - 1.0d / (4.0d * df)) - ncp) / Math.Sqrt(1.0d + q * q / (2.0d * df));
+                cdfResult = Normal.StandardCDF(Z);
+            }
+            else
+            {
+                var NCTDist = new NoncentralT(df, ncp);
+                cdfResult = NCTDist.CDF(q);
+            }
+            ANS = 1.0d - cdfResult;
             return ANS;
         }
 
diff --git a/Numerics/Distributions/Univariate/Binomial.cs b/Numerics/Distributions/Univariate/Binomial.cs
@@ -282,7 +282,7 @@ public override double CDF(double k)
             k = Math.Floor(k);
             if (k < Minimum)
                 return 0.0d;
-            if (k > Maximum)
+            if (k >= Maximum)
                 return 1.0d;
             return Beta.Incomplete(NumberOfTrials - k, k + 1d, Complement);
         }
@@ -299,7 +299,7 @@ public override double InverseCDF(double probability)
                 return Maximum;
             // Validate parameters
             if (_parametersValid == false)
-                ValidateParameters([probability, NumberOfTrials], true);
+                ValidateParameters([ProbabilityOfSuccess, NumberOfTrials], true);
             double k = 0d;
             for (int i = 0; i <= NumberOfTrials; i++)
             {
diff --git a/Numerics/Distributions/Univariate/Mixture.cs b/Numerics/Distributions/Univariate/Mixture.cs
@@ -761,13 +761,14 @@ double logLH(double[] x)
                 // Perform the expectation step
                 newLogLH = EStep(mleParameters);
 
-                // Perform the maximization step
-                mleParameters = MStep(mleParameters);
-
-                // Check convergence after M-step
+                // Check convergence before M-step to avoid pushing parameters
+                // into a degenerate state after the log-likelihood has already converged.
                 if (Math.Abs((oldLogLH - newLogLH) / oldLogLH) < Tolerance)
                     break;
 
+                // Perform the maximization step
+                mleParameters = MStep(mleParameters);
+
                 // Update log-likelihood state
                 oldLogLH = newLogLH;
 
diff --git a/Numerics/Distributions/Univariate/NoncentralT.cs b/Numerics/Distributions/Univariate/NoncentralT.cs
@@ -74,13 +74,13 @@ public NoncentralT(double degreesOfFreedom, double noncentrality)
             SetParameters(degreesOfFreedom, noncentrality);
         }
 
-        private int _degreesOfFreedom;
+        private double _degreesOfFreedom;
         private double _noncentrality;
        
         /// <summary>
         /// Gets and sets the degrees of freedom ν (nu) of the distribution.
         /// </summary>
-        public int DegreesOfFreedom
+        public double DegreesOfFreedom
         {
             get { return _degreesOfFreedom; }
             set
@@ -260,7 +260,7 @@ public override double[] MaximumOfParameters
         /// <param name="mu">The noncentrality parameter μ (mu).</param>
         public void SetParameters(double v, double mu)
         {
-            DegreesOfFreedom = (int)v;
+            DegreesOfFreedom = v;
             Noncentrality = mu;
         }
 
@@ -429,8 +429,8 @@ private double NCTDist(double t, double df, double delta)
             // 
             // Note - ITRMAX and ERRMAX may be changed to suit one's needs.
             // 
-            const int ITRMAX = 1000;
-            const double Errmax = 0.0000001d;
+            const int ITRMAX = 10000;
+            const double Errmax = 0.000000001d;
 
             // DATA ITRMAX/100.1/, ERRMAX/1.E-06/
             // 
diff --git a/Numerics/Distributions/Univariate/StudentT.cs b/Numerics/Distributions/Univariate/StudentT.cs
@@ -108,7 +108,7 @@ public StudentT(double location, double scale, double degreesOfFreedom)
 
         private double _mu;
         private double _sigma;
-        private int _degreesOfFreedom;
+        private double _degreesOfFreedom;
 
         /// <summary>
         /// Gets and sets the location parameter µ (Mu).
@@ -139,7 +139,7 @@ public double Sigma
         /// <summary>
         /// Gets and sets the degrees of freedom ν (nu) of the distribution.
         /// </summary>
-        public int DegreesOfFreedom
+        public double DegreesOfFreedom
         {
             get { return _degreesOfFreedom; }
             set
@@ -325,7 +325,7 @@ public void SetParameters(double mu, double sigma, double v)
         {
             Mu = mu;
             Sigma = sigma;
-            DegreesOfFreedom = (int)v;
+            DegreesOfFreedom = v;
         }
 
         /// <inheritdoc/>
diff --git a/Numerics/Machine Learning/Unsupervised/GaussianMixtureModel.cs b/Numerics/Machine Learning/Unsupervised/GaussianMixtureModel.cs
@@ -296,14 +296,14 @@ private void MStep()
             for (int k = 0; k < K; k++)
             {
                 double wgt = 0d;
-                for (int i = 0; i < X.NumberOfRows; i++) 
+                for (int i = 0; i < X.NumberOfRows; i++)
                     wgt += LikelihoodMatrix[i, k];
                 Weights[k] = wgt / X.NumberOfRows;
                 for (int d = 0; d < Dimension; d++)
                 {
                     // Compute centroids
                     double sum = 0;
-                    for (int i = 0; i < X.NumberOfRows; i++) 
+                    for (int i = 0; i < X.NumberOfRows; i++)
                         sum += LikelihoodMatrix[i, k] * X[i, d];
                     Means[k, d] = sum / wgt;
                     // Compute covariance
@@ -317,6 +317,14 @@ private void MStep()
                         Sigmas[k][d, j] = sum / wgt;
                     }
                 }
+
+                // Add small regularization to the diagonal to ensure the covariance
+                // matrix remains positive-definite. When a component collapses to very
+                // few points, the covariance can become singular, causing Cholesky
+                // decomposition in the E-step to fail.
+                for (int d = 0; d < Dimension; d++)
+                    MatrixRegularization.MakeSymmetricPositiveDefinite(Sigmas[k]);
+                
             }
         }
 
diff --git a/Test_Numerics/Distributions/Univariate/Test_NoncentralT.cs b/Test_Numerics/Distributions/Univariate/Test_NoncentralT.cs
@@ -114,12 +114,12 @@ public void Test_NoncentralT_InverseCDF()
         public void Test_Construction()
         {
             var t = new NoncentralT();
-            Assert.AreEqual(10,t.DegreesOfFreedom);
-            Assert.AreEqual(0, t.Noncentrality);
+            Assert.AreEqual(10d, t.DegreesOfFreedom);
+            Assert.AreEqual(0d, t.Noncentrality);
 
             var t2 = new NoncentralT(1, 1);
-            Assert.AreEqual(1, t2.DegreesOfFreedom);
-            Assert.AreEqual(1, t2.Noncentrality);
+            Assert.AreEqual(1d, t2.DegreesOfFreedom);
+            Assert.AreEqual(1d, t2.Noncentrality);
         }
 
         /// <summary>
diff --git a/Test_Numerics/Distributions/Univariate/Test_StudentT.cs b/Test_Numerics/Distributions/Univariate/Test_StudentT.cs
@@ -60,14 +60,14 @@ public class Test_StudentT
         [TestMethod()]
         public void Test_StudentT_PDF()
         {
-            var t = new StudentT(4.2d);
+            var t = new StudentT(4d);
             double pdf = t.PDF(1.4d);
             double result = 0.138377537135553d;
-            Assert.AreEqual(pdf, result, 1E-10);
-            t = new StudentT(2.5d, 0.5d, 4.2d);
+            Assert.AreEqual(result, pdf, 1E-10);
+            t = new StudentT(2.5d, 0.5d, 4d);
             pdf = t.PDF(1.4d);
             result = 0.0516476521260042d;
-            Assert.AreEqual(pdf, result, 1E-10);
+            Assert.AreEqual(result, pdf, 1E-10);
         }
 
         /// <summary>
@@ -76,14 +76,14 @@ public void Test_StudentT_PDF()
         [TestMethod()]
         public void Test_StudentT_CDF()
         {
-            var t = new StudentT(4.2d);
+            var t = new StudentT(4d);
             double cdf = t.CDF(1.4d);
             double result = 0.882949686336585d;
-            Assert.AreEqual(cdf, result, 1E-10);
-            t = new StudentT(2.5d, 0.5d, 4.2d);
+            Assert.AreEqual(result, cdf, 1E-10);
+            t = new StudentT(2.5d, 0.5d, 4d);
             cdf = t.CDF(1.4d);
             result = 0.0463263350898173d;
-            Assert.AreEqual(cdf, result, 1E-10);
+            Assert.AreEqual(result, cdf, 1E-10);
         }
 
         /// <summary>
@@ -92,16 +92,16 @@ public void Test_StudentT_CDF()
         [TestMethod()]
         public void Test_StudentT_InverseCDF()
         {
-            var t = new StudentT(4.2d);
+            var t = new StudentT(4d);
             double cdf = t.CDF(1.4d);
             double invcdf = t.InverseCDF(cdf);
             double result = 1.4d;
-            Assert.AreEqual(invcdf, result, 1E-2);
-            t = new StudentT(2.5d, 0.5d, 4.2d);
+            Assert.AreEqual(result, invcdf, 1E-2);
+            t = new StudentT(2.5d, 0.5d, 4d);
             cdf = t.CDF(1.4d);
             invcdf = t.InverseCDF(cdf);
             result = 1.4d;
-            Assert.AreEqual(invcdf, result, 1E-2);
+            Assert.AreEqual(result, invcdf, 1E-2);
         }
 
         /// <summary>
@@ -111,14 +111,14 @@ public void Test_StudentT_InverseCDF()
         public void Test_Construction()
         {
             var t = new StudentT();
-            Assert.AreEqual(0, t.Mu);
-            Assert.AreEqual(1, t.Sigma);
-            Assert.AreEqual(10, t.DegreesOfFreedom);
+            Assert.AreEqual(0d, t.Mu);
+            Assert.AreEqual(1d, t.Sigma);
+            Assert.AreEqual(10d, t.DegreesOfFreedom);
 
             var t2 = new StudentT(10, 10, 10);
-            Assert.AreEqual(10, t2.Mu);
-            Assert.AreEqual(10, t2.Sigma);
-            Assert.AreEqual(10, t2.DegreesOfFreedom);
+            Assert.AreEqual(10d, t2.Mu);
+            Assert.AreEqual(10d, t2.Sigma);
+            Assert.AreEqual(10d, t2.DegreesOfFreedom);
         }
 
         /// <summary>
@@ -148,7 +148,7 @@ public void Test_ParametersToString()
             Assert.AreEqual("Scale (σ)", t.ParametersToString[1, 0]);
             Assert.AreEqual("Degrees of Freedom (ν)", t.ParametersToString[2, 0]);
             Assert.AreEqual("0", t.ParametersToString[0, 1]);
-            Assert.AreEqual("1", t.ParametersToString[1, 1]);
+            Assert.AreEqual("1", t.ParametersToString[1,1]);
             Assert.AreEqual("10", t.ParametersToString[2, 1]);
         }
 
@@ -160,10 +160,10 @@ public void Test_Moments()
         {
             var dist = new StudentT(10, 1, 100);
             var mom = dist.CentralMoments(1E-8);
-            Assert.AreEqual(mom[0], dist.Mean, 1E-2);
-            Assert.AreEqual(mom[1], dist.StandardDeviation, 1E-2);
-            Assert.AreEqual(mom[2], dist.Skewness, 1E-2);
-            Assert.AreEqual(mom[3], dist.Kurtosis, 1E-2);
+            Assert.AreEqual(dist.Mean, mom[0], 1E-2);
+            Assert.AreEqual(dist.StandardDeviation, mom[1], 1E-2);
+            Assert.AreEqual(dist.Skewness, mom[2], 1E-2);
+            Assert.AreEqual(dist.Kurtosis, mom[3], 1E-2);
         }
 
         /// <summary>
@@ -231,7 +231,7 @@ public void Test_Skewness()
             Assert.AreEqual(0, t.Skewness);
 
             var t2 = new StudentT(1, 1, 1);
-            Assert.AreEqual(double.NaN,t2.Skewness );
+            Assert.AreEqual(double.NaN, t2.Skewness);
         }
 
         /// <summary>
@@ -241,7 +241,7 @@ public void Test_Skewness()
         public void Test_Kurtosis()
         {
             var t = new StudentT();
-            Assert.AreEqual(4, t.Kurtosis );
+            Assert.AreEqual(4, t.Kurtosis);
 
             var t2 = new StudentT(1, 1, 4);
             Assert.AreEqual(double.PositiveInfinity, t2.Kurtosis);

Original file line number	Diff line number	Diff line change
`@@ -296,14 +296,14 @@ private void MStep()`
`296`	`296`	`for (int k = 0; k < K; k++)`
`297`	`297`	`{`
`298`	`298`	`double wgt = 0d;`
`299`		`- for (int i = 0; i < X.NumberOfRows; i++)`
	`299`	`+ for (int i = 0; i < X.NumberOfRows; i++)`
`300`	`300`	`wgt += LikelihoodMatrix[i, k];`
`301`	`301`	`Weights[k] = wgt / X.NumberOfRows;`
`302`	`302`	`for (int d = 0; d < Dimension; d++)`
`303`	`303`	`{`
`304`	`304`	`// Compute centroids`
`305`	`305`	`double sum = 0;`
`306`		`- for (int i = 0; i < X.NumberOfRows; i++)`
	`306`	`+ for (int i = 0; i < X.NumberOfRows; i++)`
`307`	`307`	`sum += LikelihoodMatrix[i, k] * X[i, d];`
`308`	`308`	`Means[k, d] = sum / wgt;`
`309`	`309`	`// Compute covariance`
`@@ -317,6 +317,14 @@ private void MStep()`
`317`	`317`	`Sigmas[k][d, j] = sum / wgt;`
`318`	`318`	`}`
`319`	`319`	`}`
	`320`	`+`
	`321`	`+ // Add small regularization to the diagonal to ensure the covariance`
	`322`	`+ // matrix remains positive-definite. When a component collapses to very`
	`323`	`+ // few points, the covariance can become singular, causing Cholesky`
	`324`	`+ // decomposition in the E-step to fail.`
	`325`	`+ for (int d = 0; d < Dimension; d++)`
	`326`	`+ MatrixRegularization.MakeSymmetricPositiveDefinite(Sigmas[k]);`
	`327`	`+`
`320`	`328`	`}`
`321`	`329`	`}`
`322`	`330`