Improved Json serialization for MCMCResults and UncertaintyAnalysisResults. Improved initialization and regularization of Gaussian Mixture Model.

Author: HadenSmith

Numerics/Distributions/Univariate/Uncertainty Analysis/UncertaintyAnalysisResults.cs

@@ -163,7 +163,8 @@ public static byte[] ToByteArray(UncertaintyAnalysisResults results)
             {
                 WriteIndented = false,
                 DefaultIgnoreCondition = JsonIgnoreCondition.WhenWritingNull,
-                IncludeFields = true
+                IncludeFields = true,
+                NumberHandling = JsonNumberHandling.AllowNamedFloatingPointLiterals
             };
             // Add custom converters for unsupported types
             options.Converters.Add(new Double2DArrayConverter());
@@ -183,7 +184,8 @@ public static byte[] ToByteArray(UncertaintyAnalysisResults results)
                 var options = new JsonSerializerOptions
                 {
                     DefaultIgnoreCondition = JsonIgnoreCondition.WhenWritingNull,
-                    IncludeFields = true
+                    IncludeFields = true,
+                    NumberHandling = JsonNumberHandling.AllowNamedFloatingPointLiterals
                 };
                 // Add custom converters for unsupported types
                 options.Converters.Add(new Double2DArrayConverter());

Numerics/Machine Learning/Unsupervised/GaussianMixtureModel.cs

@@ -186,22 +186,64 @@ public GaussianMixtureModel(Matrix X, int k)
         /// <param name="kMeansPlusPlus">Determines whether to use random initialization or to use the k-Means++ method. Default is to use k-Means++.</param>
         public void Train(int seed = -1, bool kMeansPlusPlus = true)
         {
-            // 1. Initialize clusters from k-Means 
+            // 1. Initialize clusters from k-Means
             var kMeans = new KMeans(X, K);
             kMeans.Train(seed, kMeansPlusPlus);
             Means = kMeans.Means;
 
             // Give equal weight to each cluster
-            // And initialize the covariance matrix
+            // and initialize the covariance matrix from cluster data variance.
+            // Using actual data variance (instead of a tiny constant like 1e-10)
+            // prevents the first E-step from computing likelihoods from near-degenerate
+            // Gaussians, which can cause numerical overflow or underflow.
             Weights = new double[K];
             LikelihoodMatrix = new double[X.NumberOfRows, K];
             Sigmas = new Matrix[K];
             for (int k = 0; k < K; k++)
-            {              
+            {
                 Weights[k] = 1d / K;
                 Sigmas[k] = new Matrix(Dimension);
-                for (int i = 0; i < Dimension; i++)
-                    Sigmas[k][i, i] = 1E-10;
+
+                // Compute within-cluster covariance from k-means labels
+                int clusterCount = 0;
+                for (int i = 0; i < X.NumberOfRows; i++)
+                {
+                    if (kMeans.Labels[i] == k)
+                        clusterCount++;
+                }
+
+                if (clusterCount > 1)
+                {
+                    // Compute covariance from cluster members
+                    for (int d = 0; d < Dimension; d++)
+                    {
+                        for (int j = 0; j < Dimension; j++)
+                        {
+                            double sum = 0;
+                            for (int i = 0; i < X.NumberOfRows; i++)
+                            {
+                                if (kMeans.Labels[i] == k)
+                                    sum += (X[i, d] - Means[k, d]) * (X[i, j] - Means[k, j]);
+                            }
+                            Sigmas[k][d, j] = sum / clusterCount;
+                        }
+                    }
+                }
+
+                // Ensure positive-definite: floor diagonal at fraction of overall data variance
+                for (int d = 0; d < Dimension; d++)
+                {
+                    double colVar = 0;
+                    double colMean = 0;
+                    for (int i = 0; i < X.NumberOfRows; i++)
+                        colMean += X[i, d];
+                    colMean /= X.NumberOfRows;
+                    for (int i = 0; i < X.NumberOfRows; i++)
+                        colVar += (X[i, d] - colMean) * (X[i, d] - colMean);
+                    colVar /= X.NumberOfRows;
+
+                    Sigmas[k][d, d] = Math.Max(Sigmas[k][d, d], 1E-6 * colVar);
+                }
             }
 
             // 2. Optimize clusters
@@ -318,13 +360,26 @@ private void MStep()
                     }
                 }
 
-                // 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.
+                // Floor diagonal at a fraction of the overall data variance to prevent
+                // component collapse. When a component captures very few points, its
+                // covariance can become singular, causing Cholesky decomposition in the
+                // E-step to fail. This mirrors sklearn's reg_covar parameter.
                 for (int d = 0; d < Dimension; d++)
-                    MatrixRegularization.MakeSymmetricPositiveDefinite(Sigmas[k]);
-                
+                {
+                    double colVar = 0;
+                    double colMean = 0;
+                    for (int i = 0; i < X.NumberOfRows; i++)
+                        colMean += X[i, d];
+                    colMean /= X.NumberOfRows;
+                    for (int i = 0; i < X.NumberOfRows; i++)
+                        colVar += (X[i, d] - colMean) * (X[i, d] - colMean);
+                    colVar /= X.NumberOfRows;
+
+                    Sigmas[k][d, d] = Math.Max(Sigmas[k][d, d], 1E-6 * colVar);
+                }
+
+                // Ensure the full covariance matrix remains symmetric positive-definite
+                MatrixRegularization.MakeSymmetricPositiveDefinite(Sigmas[k]);
             }
         }
 

Numerics/Sampling/MCMC/Support/MCMCResults.cs

@@ -192,7 +192,8 @@ public static byte[] ToByteArray(MCMCResults mcmcResults)
             {
                 WriteIndented = false,
                 DefaultIgnoreCondition = JsonIgnoreCondition.WhenWritingNull,
-                IncludeFields = true
+                IncludeFields = true,
+                NumberHandling = JsonNumberHandling.AllowNamedFloatingPointLiterals
             };
             options.Converters.Add(new Double2DArrayConverter());
             options.Converters.Add(new HistogramConverter());
@@ -208,7 +209,8 @@ public static byte[] ToByteArray(MCMCResults mcmcResults)
             var options = new JsonSerializerOptions
             {
                 DefaultIgnoreCondition = JsonIgnoreCondition.WhenWritingNull,
-                IncludeFields = true
+                IncludeFields = true,
+                NumberHandling = JsonNumberHandling.AllowNamedFloatingPointLiterals
             };
             options.Converters.Add(new Double2DArrayConverter());
             options.Converters.Add(new HistogramConverter());