Multiple edge case bug fixes. Added Yeo-Johnson transformation, Unimodality test, new methods in TimeSeries class, and fixed minor bugs in distributions and MCMC sampling methods.

Author: HadenSmith

4.8.1/Numerics/Data/Statistics/BoxCox.cs

@@ -52,31 +52,23 @@ public class BoxCox
         /// Fit the transformation parameters using maximum likelihood estimation.
         /// </summary>
         /// <param name="values">The list of values to transform.</param>
-        /// <param name="lambda1">Output. The transformation exponent. Range -5 to +5.</param>
-        /// <param name="lambda2">Output. The transformation shift for negative values.</param>
+        /// <param name="lambda">Output. The transformation exponent. Range -5 to +5.</param>
         /// <remarks>
         /// https://www.rdocumentation.org/packages/EnvStats/versions/2.4.0/topics/boxcox
         /// </remarks>
-        public static void FitLambda(IList<double> values, out double lambda1, out double lambda2)
+        public static void FitLambda(IList<double> values, out double lambda)
         {
-            int n = values.Count;
-            double l1 = 0d;
-            double l2 = 0d;
-            double min = Statistics.Minimum(values);
-            if (min <= 0d) l2 = 0d - min + Tools.DoubleMachineEpsilon;
-            
-            Func<double, double> func = lambda =>
-            {
-                return LogLikelihood(values, lambda, l2);
-            };
-
             // Solve with Brent 
-            var brent = new BrentSearch(func, -5d, 5d);
+            var brent = new BrentSearch((x) => { return LogLikelihood(values, x); }, -5d, 5d);
             brent.Maximize();
-            l1 = brent.BestParameterSet.Values[0];
-            // Set parameters
-            lambda1 = l1;
-            lambda2 = l2;
+            if (brent.Status == OptimizationStatus.Success)
+            {
+                lambda = brent.BestParameterSet.Values[0];
+            }
+            else
+            {
+                lambda = double.NaN;
+            }
         }
 
         /// <summary>
@@ -85,21 +77,20 @@ public static void FitLambda(IList<double> values, out double lambda1, out doubl
         /// </summary>
         /// <param name="values">The list of values to transform.</param>
         /// <param name="lambda1">The transformation exponent. Range -5 to +5.</param>
-        /// <param name="lambda2">The transformation shift for negative values.</param>
         /// <returns>
         /// The value of log-likelihood function evaluated at the given values and lambdas.
         /// </returns>
-        public static double LogLikelihood(IList<double> values, double lambda1, double lambda2)
+        public static double LogLikelihood(IList<double> values, double lambda1)
         {
             int n = values.Count;
             var y = new double[n];
             double mu = 0d;
             var sumX = 0d;
             for (int i = 0; i < n; i++)
             {
-                y[i] = Transform(values[i], lambda1, lambda2);
+                y[i] = Transform(values[i], lambda1);
                 mu += y[i];
-                sumX += Math.Log(values[i] + lambda2);
+                sumX += Math.Log(values[i]);
             }
             mu = mu / n;
             double sse = 0d;
@@ -110,57 +101,78 @@ public static double LogLikelihood(IList<double> values, double lambda1, double
             return ll;
         }
 
+        /// <summary>
+        /// Computes the Log-Jacobian used to adjust the log-likelihood function. 
+        /// </summary>
+        /// <param name="values">The list of values to transform.</param>
+        /// <param name="lambda">The transformation exponent. Range -5 to +5.</param>
+        /// <returns>Returns the Log Jacobian.</returns>
+        public static double LogJacobian(IList<double> values, double lambda)
+        {
+            double logJacobianSum = 0d;
+            int n = values.Count;
+
+            for (int i = 0; i < n; i++)
+            {
+                double xi = values[i];
+
+                if (xi <= 0)
+                    return double.NegativeInfinity; // Box-Cox undefined for non-positive values
+
+                // For Box-Cox: log derivative is (lambda - 1) * log(x)
+                logJacobianSum += (lambda - 1d) * Math.Log(xi);
+            }
+
+            return logJacobianSum;
+        }
+
         /// <summary>
         /// Returns the Box-Cox transformation of the value.
         /// </summary>
         /// <param name="value">The value to transform.</param>
-        /// <param name="lambda1">The transformation exponent. Range -5 to +5.</param>
-        /// <param name="lambda2">The transformation shift for negative values.</param>
-        public static double Transform(double value, double lambda1, double lambda2 = 0.0d)
+        /// <param name="lambda">The transformation exponent. Range -5 to +5.</param>
+        public static double Transform(double value, double lambda)
         {
-            if (Math.Abs(lambda1) > 5d) return double.NaN;
-            if (lambda1 == 0d) return Math.Log(value + lambda2);
-            return (Math.Pow(value + lambda2, lambda1) - 1.0d) / lambda1;
+            if (Math.Abs(lambda) > 5d) return double.NaN;
+            if (Math.Abs(lambda) < 1e-8) return Math.Log(value);
+            return (Math.Pow(value, lambda) - 1.0d) / lambda;
         }
 
         /// <summary>
         /// Returns the Box-Cox transformation of each value in the list.
         /// </summary>
         /// <param name="values">The list of values to transform.</param>
-        /// <param name="lambda1">The transformation exponent. Range -5 to +5.</param>
-        /// <param name="lambda2">The transformation shift for negative values.</param>
-        public static List<double> Transform(IList<double> values, double lambda1, double lambda2 = 0.0d)
+        /// <param name="lambda">The transformation exponent. Range -5 to +5.</param>
+        public static List<double> Transform(IList<double> values, double lambda)
         {
             var newValues = new List<double>();
             for (int i = 0; i < values.Count; i++)
-                newValues.Add(Transform(values[i], lambda1, lambda2));
+                newValues.Add(Transform(values[i], lambda));
             return newValues;
         }
 
         /// <summary>
         /// Returns the reverse of the Box-Cox transformed value.
         /// </summary>
         /// <param name="value">The value to reverse transform.</param>
-        /// <param name="lambda1">The transformation exponent. Range -5 to +5.</param>
-        /// <param name="lambda2">The transformation shift for negative values.</param>
-        public static double ReverseTransform(double value, double lambda1, double lambda2 = 0.0d)
+        /// <param name="lambda">The transformation exponent. Range -5 to +5.</param>
+        public static double InverseTransform(double value, double lambda)
         {
-            if (Math.Abs(lambda1) > 5d) return double.NaN;
-            if (lambda1 == 0d) return Math.Exp(value) - lambda2;
-            return Math.Pow(value * lambda1 + 1.0d, 1.0d / lambda1) - lambda2;
+            if (Math.Abs(lambda) > 5d) return double.NaN;
+            if (Math.Abs(lambda) < 1e-8) return Math.Exp(value);
+            return Math.Pow(value * lambda + 1.0d, 1.0d / lambda);
         }
 
         /// <summary>
-        /// Returns the reverse of each Box-Cox transformed value in the list.
+        /// Returns the inverse of each Box-Cox transformed value in the list.
         /// </summary>
         /// <param name="values">The list of values to reverse transform.</param>
-        /// <param name="lambda1">The transformation exponent. Range -5 to +5.</param>
-        /// <param name="lambda2">The transformation shift for negative values.</param>
-        public static List<double> ReverseTransform(IList<double> values, double lambda1, double lambda2 = 0.0d)
+        /// <param name="lambda">The transformation exponent. Range -5 to +5.</param>
+        public static List<double> InverseTransform(IList<double> values, double lambda)
         {
             var newValues = new List<double>();
             for (int i = 0; i < values.Count; i++)
-                newValues.Add(ReverseTransform(values[i], lambda1, lambda2));
+                newValues.Add(InverseTransform(values[i], lambda));
             return newValues;
         }
     }

4.8.1/Numerics/Data/Statistics/Correlation.cs

@@ -117,8 +117,8 @@ public static double Spearman(IList<double> sample1, IList<double> sample2)
                 throw new ArgumentOutOfRangeException(nameof(sample1), "The sample arrays must be the same length.");
             }
 
-            var rank1 = Statistics.RanksInplace(sample1.ToArray());
-            var rank2 = Statistics.RanksInplace(sample2.ToArray());
+            var rank1 = Statistics.RanksInPlace(sample1.ToArray());
+            var rank2 = Statistics.RanksInPlace(sample2.ToArray());
             return Pearson(rank1, rank2);
         }
 

4.8.1/Numerics/Data/Statistics/GoodnessOfFit.cs

@@ -113,28 +113,6 @@ public static double[] AICWeights(IList<double> aicValues)
             return weights;
         }
 
-        /// <summary>
-        /// Returns an array of weights based on a list of model BIC values.
-        /// </summary>
-        /// <param name="bicValues">The list of model BIC values.</param>
-        public static double[] BICWeights(IList<double> bicValues)
-        {
-            var min = Tools.Min(bicValues);
-            var weights = new double[bicValues.Count];
-            var num = new double[bicValues.Count];
-            double sum = 0;
-            for (int i = 0; i < bicValues.Count; i++)
-            {
-                num[i] = Math.Exp(-0.5 * (bicValues[i] - min));
-                sum += num[i];
-            }
-            for (int i = 0; i < bicValues.Count; i++)
-            {
-                weights[i] = num[i] / sum;
-            }
-            return weights;
-        }
-
         /// <summary>
         /// Gets the Root Mean Square Error (RMSE) of the model compared to the observed data.
         /// </summary>

4.8.1/Numerics/Data/Statistics/HypothesisTests.cs

@@ -32,6 +32,7 @@
 using System.Collections.Generic;
 using System.Linq;
 using Numerics.Distributions;
+using Numerics.MachineLearning;
 using Numerics.Mathematics.LinearAlgebra;
 using Numerics.Mathematics.SpecialFunctions;
 
@@ -142,7 +143,7 @@ public static double PairedTtest(IList<double> sample1, IList<double> sample2)
         /// <param name="sample2">Data sample 2.</param>
         /// <returns>Returns the p-value of the test statistic.</returns>
         public static double Ftest(IList<double> sample1, IList<double> sample2)
-        {        
+        {
             var meanVar1 = Statistics.MeanVariance(sample1);
             var meanVar2 = Statistics.MeanVariance(sample2);
             int n1 = sample1.Count, n2 = sample2.Count;
@@ -224,7 +225,7 @@ public static double WaldWolfowitzTest(IList<double> sample)
             R += xn * x1;
             double s12 = Tools.Sqr(S1);
             double s22 = Tools.Sqr(S2);
-            double s14 = Tools.Pow(S1,4);
+            double s14 = Tools.Pow(S1, 4);
 
             eR = (s12 - S2) / (n - 1);
             varR = (s22 - S4) / (n - 1) - Tools.Sqr(eR);
@@ -294,7 +295,7 @@ public static double MannKendallTest(IList<double> sample)
             int i, j, n = sample.Count;
             if (n < 10) throw new ArgumentException("The sample size must be greater than or equal to 10.");
 
-            double S = 0, T = 0, varS, z; 
+            double S = 0, T = 0, varS, z;
 
             for (i = 0; i < n - 1; i++)
             {
@@ -331,5 +332,40 @@ public static double LinearTrendTest(IList<double> indices, IList<double> sample
             return (1 - tdist.CDF(Math.Abs(lm.Parameters[1] / lm.ParameterStandardErrors[1]))) * 2;
         }
 
+        /// <summary>
+        /// The Gaussian Mixture Model Unimodality Test.
+        /// </summary>
+        /// <param name="sample">Time series sample data.</param>
+        /// <returns>Returns the p-value of the test statistic.</returns>
+        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();
+            var logLH1 = gmm1.LogLikelihood;
+
+            // Fit 2-component GMM (potentially bimodal)
+            var gmm2 = new GaussianMixtureModel(sample.ToArray(), 2);
+            gmm2.Train();
+            var logLH2 = gmm2.LogLikelihood;
+
+            // 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);
+
+            return pval;
+
+        }
+
+
     }
+
+
 }

4.8.1/Numerics/Data/Statistics/MultipleGrubbsBeckTest.cs

@@ -199,8 +199,7 @@ private static double GGBCRITP(int N, int R, double ETA)
             // The number of low outliers computed by this method is consistent with the results from the FORTRAN code.
             // Further testing is necessary to identify any edge cases where the adaptive Simpson's rule might prove insufficient.
             var sr = new AdaptiveSimpsonsRule(FGGB, 1E-16, 1 - 1E-16);
-            sr.RelativeTolerance = 1E-4;
-            sr.MaxDepth = 50;
+            sr.MaxDepth = 25;
             sr.ReportFailure = false;
             sr.Integrate();
             return sr.Status != IntegrationStatus.Failure ? sr.Result : double.NaN;