Bug fixes and enhancements

.github/workflows/Snapshot.yml

@@ -11,7 +11,8 @@ jobs:
       dotnet-version: '10.0.x'
       project-names: 'Numerics'
       version: '2.0.0'
-      run-tests: true
+      # PR integration already runs the full multi-target test suite.
+      run-tests: false
       nuget-source: 'https://www.hec.usace.army.mil/nexus/repository/consequences-nuget-public/'
     secrets:
       NUGET_API_KEY: ${{ secrets.NEXUS_NUGET_APIKEY }}

Numerics/Data/Statistics/YeoJohnson.cs

@@ -37,6 +37,20 @@ public static void FitLambda(IList<double> values, out double lambda)
             }
         }
 
+        /// <summary>
+        /// Fit the transform parameter using maximum likelihood estimation.
+        /// </summary>
+        /// <param name="values">The list of values to transform.</param>
+        /// <returns>The fitted transformation exponent.</returns>
+        public static double FitLambda(IList<double> values)
+        {
+            if (values == null) throw new ArgumentNullException(nameof(values));
+            if (values.Count < 2) throw new ArgumentException("At least 2 values are required to fit lambda.", nameof(values));
+
+            FitLambda(values, out double lambda);
+            return lambda;
+        }
+
         /// <summary>
         /// The log-likelihood function. The transformed observations are assumed to come from a
         /// normal distribution. The change of variable formula is used to write the log-likelihood function.
@@ -160,6 +174,19 @@ public static double Transform(double value, double lambda)
             }
         }
 
+        /// <summary>
+        /// Returns the derivative of the Yeo-Johnson transformation with respect to the original value.
+        /// </summary>
+        /// <param name="value">The value at which to evaluate the derivative.</param>
+        /// <param name="lambda">The transformation exponent. Range -5 to +5.</param>
+        public static double Derivative(double value, double lambda)
+        {
+            if (Math.Abs(lambda) > 5d) return double.NaN;
+            return value >= 0d
+                ? Math.Pow(value + 1d, lambda - 1d)
+                : Math.Pow(-value + 1d, 1d - lambda);
+        }
+
         /// <summary>
         /// Returns the Yeo-Johnson transformation of each value in the list.
         /// </summary>

Numerics/Distributions/Bivariate Copulas/Base/BivariateCopula.cs

@@ -110,6 +110,8 @@ public bool ParametersValid
         public double LogPDF(double u, double v)
         {
             double f = PDF(u, v);
+            if (double.IsNaN(f) || double.IsInfinity(f) || f <= 0d)
+                return double.NegativeInfinity;
             return Math.Log(f);
         }
 
@@ -194,7 +196,7 @@ public double PseudoLogLikelihood(IList<double> sampleDataX, IList<double> sampl
             double LogLH = 0;
             for (int i = 0; i < sampleDataX.Count; i++)
                 LogLH += LogPDF(sampleDataX[i], sampleDataY[i]);
-            if (double.IsNaN(LogLH) || double.IsInfinity(LogLH)) return double.MinValue;
+            if (double.IsNaN(LogLH) || double.IsInfinity(LogLH)) return double.NegativeInfinity;
             return LogLH;
         }
 
@@ -212,7 +214,7 @@ public double IFMLogLikelihood(IList<double> sampleDataX, IList<double> sampleDa
             double LogLH = 0;
             for (int i = 0; i < sampleDataX.Count; i++)
                 LogLH += LogPDF(MarginalDistributionX.CDF(sampleDataX[i]), MarginalDistributionY.CDF(sampleDataY[i]));
-            if (double.IsNaN(LogLH) || double.IsInfinity(LogLH)) return double.MinValue;
+            if (double.IsNaN(LogLH) || double.IsInfinity(LogLH)) return double.NegativeInfinity;
             return LogLH;
         }
 
@@ -230,7 +232,7 @@ public double LogLikelihood(IList<double> sampleDataX, IList<double> sampleDataY
             double LogLH = 0;
             for (int i = 0; i < sampleDataX.Count; i++)
                 LogLH += LogPDF(MarginalDistributionX.CDF(sampleDataX[i]), MarginalDistributionY.CDF(sampleDataY[i])) + MarginalDistributionX.LogPDF(sampleDataX[i]) + MarginalDistributionY.LogPDF(sampleDataY[i]);
-            if (double.IsNaN(LogLH) || double.IsInfinity(LogLH)) return double.MinValue;
+            if (double.IsNaN(LogLH) || double.IsInfinity(LogLH)) return double.NegativeInfinity;
             return LogLH;
         }
 

Numerics/Distributions/Multivariate/Base/MultivariateDistribution.cs

@@ -55,8 +55,8 @@ public abstract class MultivariateDistribution : IMultivariateDistribution
         public virtual double LogPDF(double[] x)
         {
             double f = PDF(x);
-            // If the PDF returns an invalid probability, then return the worst log-probability.
-            if (double.IsNaN(f) || double.IsInfinity(f) || f <= 0d) return double.MinValue;
+            // If the PDF returns an invalid probability, then return log(0).
+            if (double.IsNaN(f) || double.IsInfinity(f) || f <= 0d) return double.NegativeInfinity;
             return Math.Log(f);
         }
 
@@ -73,8 +73,8 @@ public virtual double LogPDF(double[] x)
         public virtual double LogCDF(double[] x)
         {
             double F = CDF(x);
-            // If the CDF returns an invalid probability, then return the worst log-probability.
-            if (double.IsNaN(F) || double.IsInfinity(F) || F <= 0d) return double.MinValue;
+            // If the CDF returns an invalid probability, then return log(0).
+            if (double.IsNaN(F) || double.IsInfinity(F) || F <= 0d) return double.NegativeInfinity;
             return Math.Log(F);
         }
 
@@ -94,9 +94,9 @@ public virtual double CCDF(double[] x)
         public virtual double LogCCDF(double[] x)
         {
             double cF = CCDF(x);
-            // If the CCDF returns an invalid probability, then return the worst log-probability.
-            if (double.IsNaN(cF) || cF <= 0d)
-                return int.MinValue;
+            // If the CCDF returns an invalid probability, then return log(0).
+            if (double.IsNaN(cF) || double.IsInfinity(cF) || cF <= 0d)
+                return double.NegativeInfinity;
             return Math.Log(cF);
         }
 
@@ -106,4 +106,4 @@ public virtual double LogCCDF(double[] x)
         public abstract MultivariateDistribution Clone();
 
     }
-}
+}

Numerics/Distributions/Multivariate/Dirichlet.cs

@@ -282,7 +282,7 @@ public override double PDF(double[] x)
         /// </summary>
         /// <param name="x">A point on the simplex. Components must be positive and sum to 1.</param>
         /// <returns>
-        /// The log-density at x. Returns <see cref="double.MinValue"/> if x is not on the simplex.
+        /// The log-density at x. Returns <see cref="double.NegativeInfinity"/> if x is not on the simplex.
         /// </returns>
         /// <remarks>
         /// <para>
@@ -295,16 +295,16 @@ public override double PDF(double[] x)
         /// </remarks>
         public override double LogPDF(double[] x)
         {
-            if (x == null || x.Length != _dimension) return double.MinValue;
+            if (x == null || x.Length != _dimension) return double.NegativeInfinity;
 
             // Check that x is on the simplex
             double sum = 0;
             for (int i = 0; i < _dimension; i++)
             {
-                if (x[i] <= 0 || x[i] > 1) return double.MinValue;
+                if (x[i] <= 0 || x[i] > 1) return double.NegativeInfinity;
                 sum += x[i];
             }
-            if (Math.Abs(sum - 1.0) > 1e-10) return double.MinValue;
+            if (Math.Abs(sum - 1.0) > 1e-10) return double.NegativeInfinity;
 
             double logPdf = -_logNormalization;
             for (int i = 0; i < _dimension; i++)

Numerics/Distributions/Multivariate/Multinomial.cs

@@ -204,20 +204,20 @@ public override double PDF(double[] x)
         /// Computes the log of the probability mass function.
         /// </summary>
         /// <param name="x">The count vector. Each xᵢ must be a non-negative integer and Σxᵢ = N.</param>
-        /// <returns>The log probability. Returns <see cref="double.MinValue"/> if x is not valid.</returns>
+        /// <returns>The log probability. Returns <see cref="double.NegativeInfinity"/> if x is not valid.</returns>
         public double LogPMF(double[] x)
         {
-            if (x == null || x.Length != _dimension) return double.MinValue;
+            if (x == null || x.Length != _dimension) return double.NegativeInfinity;
 
             // Validate: all non-negative integers that sum to N
             int sum = 0;
             for (int i = 0; i < _dimension; i++)
             {
                 int xi = (int)Math.Round(x[i]);
-                if (xi < 0 || Math.Abs(x[i] - xi) > 1e-10) return double.MinValue;
+                if (xi < 0 || Math.Abs(x[i] - xi) > 1e-10) return double.NegativeInfinity;
                 sum += xi;
             }
-            if (sum != _n) return double.MinValue;
+            if (sum != _n) return double.NegativeInfinity;
 
             // Compute in log-space: log(N!) - sum(log(xi!)) + sum(xi * log(pi))
             double logPmf = Factorial.LogFactorial(_n);
@@ -227,7 +227,7 @@ public double LogPMF(double[] x)
                 logPmf -= Factorial.LogFactorial(xi);
                 if (xi > 0)
                 {
-                    if (_p[i] <= 0) return double.MinValue;
+                    if (_p[i] <= 0) return double.NegativeInfinity;
                     logPmf += xi * Math.Log(_p[i]);
                 }
             }

Numerics/Distributions/Multivariate/MultivariateNormal.cs

@@ -338,7 +338,7 @@ public override double PDF(double[] x)
         public override double LogPDF(double[] x)
         {
             double f = -0.5d * Mahalanobis(x) + _lnconstant;
-            if (double.IsNaN(f) || double.IsInfinity(f)) return double.MinValue;
+            if (double.IsNaN(f) || double.IsInfinity(f)) return double.NegativeInfinity;
             return f;
         }
 
@@ -2149,4 +2149,4 @@ public override MultivariateDistribution Clone()
         }
 
     }
-}
+}

Numerics/Distributions/Multivariate/MultivariateStudentT.cs

@@ -400,7 +400,7 @@ public override double LogPDF(double[] x)
         {
             double Q = Mahalanobis(x);
             double f = _lnconstant - ((_degreesOfFreedom + _dimension) / 2.0) * Math.Log(1.0 + Q / _degreesOfFreedom);
-            if (double.IsNaN(f) || double.IsInfinity(f)) return double.MinValue;
+            if (double.IsNaN(f) || double.IsInfinity(f)) return double.NegativeInfinity;
             return f;
         }
 
@@ -705,4 +705,4 @@ public override MultivariateDistribution Clone()
         }
 
     }
-}
+}

Numerics/Distributions/Univariate/Base/UnivariateDistributionBase.cs

@@ -150,7 +150,7 @@ public double LogLikelihood(IList<double> sample)
             {        
                 logLH += LogPDF(sample[i]);
             }
-            if (double.IsNaN(logLH) || double.IsInfinity(logLH)) return double.MinValue;
+            if (double.IsNaN(logLH) || double.IsInfinity(logLH)) return double.NegativeInfinity;
             return logLH;
         }
 
@@ -204,6 +204,8 @@ public double LogLikelihood_Intervals(double lowerLimit, double upperLimit)
         public virtual double LogPDF(double x)
         {
             double f = PDF(x);
+            if (double.IsNaN(f) || double.IsInfinity(f) || f <= 0.0)
+                return double.NegativeInfinity;
             return Math.Log(f);
         }
 
@@ -220,6 +222,8 @@ public virtual double HF(double x)
         public virtual double LogCDF(double x)
         {
             double F = CDF(x);
+            if (double.IsNaN(F) || double.IsInfinity(F) || F <= 0.0)
+                return double.NegativeInfinity;
             return Math.Log(F);
         }
 
@@ -233,6 +237,8 @@ public virtual double CCDF(double x)
         public virtual double LogCCDF(double x)
         {
             double cF = CCDF(x);
+            if (double.IsNaN(cF) || double.IsInfinity(cF) || cF <= 0.0)
+                return double.NegativeInfinity;
             return Math.Log(cF);
         }
 

Numerics/Distributions/Univariate/CompetingRisks.cs

@@ -58,8 +58,8 @@ public CompetingRisks(IUnivariateDistribution[] distributions)
         private bool _mvnCreated = false;
         private MultivariateNormal _mvn = null!;
 
-        // Minimum log value to prevent -Infinity in log-likelihood
-        private const double _logZero = -745.0; // ≈ log(double.MinValue that's positive)
+        // Soft finite floor used in tail arithmetic before returning the final log-density.
+        private const double _logZero = -745.0;
         private const double _minDensity = 1E-300;
 
         /// <summary>
@@ -579,7 +579,7 @@ public override double LogPDF(double x)
                 // For dependent cases, fall back to numerical differentiation
                 // but use a more stable approach
                 double pdf = NumericalDerivative.Derivative(CDF, x);
-                return pdf > _minDensity ? Math.Log(pdf) : double.MinValue;
+                return pdf > _minDensity ? Math.Log(pdf) : double.NegativeInfinity;
             }
 
         }
@@ -648,7 +648,7 @@ private double LogPDFMinimumIndependent(double x)
             // Final result: log(f) = log(sum h_i) + sum log(S_i)
             double logPdf = logSumHazard + sumLogSurvival;
 
-            return double.IsNaN(logPdf) || double.IsInfinity(logPdf) ? double.MinValue : logPdf;
+            return double.IsNaN(logPdf) || double.IsInfinity(logPdf) ? double.NegativeInfinity : logPdf;
         }
 
         /// <summary>
@@ -714,7 +714,7 @@ private double LogPDFMaximumIndependent(double x)
             // Final result: log(f) = log(sum f_i/F_i) + sum log(F_i)
             double logPdf = logSumRatio + sumLogCdf;
 
-            return double.IsNaN(logPdf) || double.IsInfinity(logPdf) ? double.MinValue : logPdf;
+            return double.IsNaN(logPdf) || double.IsInfinity(logPdf) ? double.NegativeInfinity : logPdf;
         }
 
 

Numerics/Distributions/Univariate/Mixture.cs

@@ -672,14 +672,21 @@ double EStep(double[] x)
                 for (int i = 0; i < N; i++)
                 {
                     // Get max likelihood
-                    double max = double.MinValue;
+                    double max = double.NegativeInfinity;
                     for (int k = 0; k < K; k++)
                     {
                         if (likelihood[i, k] > max)
                         {
                             max = likelihood[i, k];
                         }
                     }
+                    if (double.IsNegativeInfinity(max))
+                    {
+                        for (int k = 0; k < K; k++)
+                            likelihood[i, k] = 0.0;
+                        return double.NegativeInfinity;
+                    }
+
                     // log-sum-exp trick begins here
                     double sum = 0;
                     for (int k = 0; k < K; k++)
@@ -721,7 +728,7 @@ double logLH(double[] x)
                 var dist = (Mixture)Clone();
                 dist.SetParameters(mleWeights, x);
                 double lh = dist.LogLikelihood(sample);
-                if (double.IsNaN(lh) || double.IsInfinity(lh)) return double.MinValue;
+                if (double.IsNaN(lh) || double.IsInfinity(lh)) return double.NegativeInfinity;
                 return lh;
             }
 

Numerics/Functions/Link Functions/FisherZLink.cs

@@ -0,0 +1,61 @@
+using System;
+using System.Xml.Linq;
+
+namespace Numerics.Functions
+{
+    /// <summary>
+    /// Fisher z link function mapping correlations from (-1, 1) to the unconstrained real line.
+    /// </summary>
+    /// <remarks>
+    /// <para>
+    ///     Domain: -1 &lt; x &lt; 1. h(x) = atanh(x), h&#8315;&#185;(&#951;) = tanh(&#951;),
+    ///     h&#8242;(x) = 1 / (1 - x^2).
+    /// </para>
+    /// <para>
+    ///     This link is commonly used for correlation parameters and other signed bounded parameters.
+    /// </para>
+    /// <para>
+    ///     <b> Authors: </b>
+    ///     Haden Smith, USACE Risk Management Center, cole.h.smith@usace.army.mil
+    /// </para>
+    /// </remarks>
+    [Serializable]
+    public sealed class FisherZLink : ILinkFunction
+    {
+        /// <inheritdoc/>
+        /// <exception cref="ArgumentOutOfRangeException">Thrown when <paramref name="x"/> is not in (-1, 1).</exception>
+        public double Link(double x)
+        {
+            if (!Tools.IsFinite(x) || x <= -1d || x >= 1d)
+                throw new ArgumentOutOfRangeException(nameof(x), "FisherZLink requires -1 < x < 1.");
+
+            return 0.5d * Math.Log((1d + x) / (1d - x));
+        }
+
+        /// <inheritdoc/>
+        public double InverseLink(double eta)
+        {
+            if (eta >= 0d)
+            {
+                double e = Math.Exp(-2d * eta);
+                return (1d - e) / (1d + e);
+            }
+
+            double expEta = Math.Exp(2d * eta);
+            return (expEta - 1d) / (expEta + 1d);
+        }
+
+        /// <inheritdoc/>
+        /// <exception cref="ArgumentOutOfRangeException">Thrown when <paramref name="x"/> is not in (-1, 1).</exception>
+        public double DLink(double x)
+        {
+            if (!Tools.IsFinite(x) || x <= -1d || x >= 1d)
+                throw new ArgumentOutOfRangeException(nameof(x), "FisherZLink derivative requires -1 < x < 1.");
+
+            return 1d / (1d - x * x);
+        }
+
+        /// <inheritdoc/>
+        public XElement ToXElement() => new XElement(nameof(FisherZLink));
+    }
+}