diff --git a/.github/workflows/Snapshot.yml b/.github/workflows/Snapshot.yml index ea5276b1..6a7bb545 100644 --- a/.github/workflows/Snapshot.yml +++ b/.github/workflows/Snapshot.yml @@ -10,8 +10,8 @@ jobs: with: dotnet-version: '10.0.x' project-names: 'Numerics' - # Snapshot workflow appends .-dev, so this tracks the next development patch after v2.1.2. - version: '2.1.3' + # Snapshot workflow appends .-dev, so this tracks the next development patch after v2.1.3. + version: '2.1.4' # 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/' diff --git a/.gitignore b/.gitignore index 0049b80c..69cc224c 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,4 @@ -## Ignore Visual Studio temporary files, build results, and +## Ignore Visual Studio temporary files, build results, and ## files generated by popular Visual Studio add-ons. ## ## Get latest from https://github.com/github/gitignore/blob/master/VisualStudio.gitignore @@ -370,4 +370,11 @@ MigrationBackup/ #Numerics Specific /TestResults -/.claude/settings.local.json + +# Local developer guidance and settings. These files are for developer machines only. +/CLAUDE.md +/Claude.md +/AGENTS.md +/Agents.md +/.claude/ +/.agents/ diff --git a/CITATION.cff b/CITATION.cff index c3bf0a23..7ae74765 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -2,8 +2,8 @@ cff-version: 1.2.0 message: "If you use this software, please cite our article in the Journal of Open Source Software." type: software title: "Numerics: A .NET Library for Numerical Computing, Statistical Analysis, and Risk Assessment" -version: "2.1.2" -date-released: "2026-07-09" +version: "2.1.3" +date-released: "2026-07-15" license: 0BSD repository-code: "https://github.com/USACE-RMC/Numerics" url: "https://github.com/USACE-RMC/Numerics" diff --git a/Directory.Build.props b/Directory.Build.props new file mode 100644 index 00000000..ece0017c --- /dev/null +++ b/Directory.Build.props @@ -0,0 +1,13 @@ + + + + $(MSBuildThisFileDirectory) + true + + + + true + $(WarningsAsErrors);CS1570;CS1571;CS1572;CS1573;CS1574;CS1584;CS1587;CS1589;CS1591 + + + diff --git a/Numerics/Data/Interpolation/Bilinear.cs b/Numerics/Data/Interpolation/Bilinear.cs index 34d8ed3f..096f18ea 100644 --- a/Numerics/Data/Interpolation/Bilinear.cs +++ b/Numerics/Data/Interpolation/Bilinear.cs @@ -142,11 +142,11 @@ public double Interpolate(double x1, double x2) x1ub = X1Values[X1LI.Count - 1]; if (X1Transform == Transform.Logarithmic) { - x1 = Math.Log10(x1); - x1i = Math.Log10(x1i); - x1ii = Math.Log10(x1ii); - x1lb = Math.Log10(x1lb); - x1ub = Math.Log10(x1ub); + x1 = Tools.Log10(x1); + x1i = Tools.Log10(x1i); + x1ii = Tools.Log10(x1ii); + x1lb = Tools.Log10(x1lb); + x1ub = Tools.Log10(x1ub); } else if (X1Transform == Transform.NormalZ) { @@ -164,11 +164,11 @@ public double Interpolate(double x1, double x2) x2ub = X2Values[X2LI.Count - 1]; if (X2Transform == Transform.Logarithmic) { - x2 = Math.Log10(x2); - x2j = Math.Log10(x2j); - x2jj = Math.Log10(x2jj); - x2lb = Math.Log10(x2lb); - x2ub = Math.Log10(x2ub); + x2 = Tools.Log10(x2); + x2j = Tools.Log10(x2j); + x2jj = Tools.Log10(x2jj); + x2lb = Tools.Log10(x2lb); + x2ub = Tools.Log10(x2ub); } else if (X2Transform == Transform.NormalZ) { @@ -186,10 +186,10 @@ public double Interpolate(double x1, double x2) yijj = YValues[i, j + 1]; if (YTransform == Transform.Logarithmic) { - yij = Math.Log10(yij); - yiij = Math.Log10(yiij); - yiijj = Math.Log10(yiijj); - yijj = Math.Log10(yijj); + yij = Tools.Log10(yij); + yiij = Tools.Log10(yiij); + yiijj = Tools.Log10(yiijj); + yijj = Tools.Log10(yijj); } else if (YTransform == Transform.NormalZ) { diff --git a/Numerics/Data/Interpolation/Support/Search.cs b/Numerics/Data/Interpolation/Support/Search.cs index 9c5f43e6..21cee576 100644 --- a/Numerics/Data/Interpolation/Support/Search.cs +++ b/Numerics/Data/Interpolation/Support/Search.cs @@ -398,16 +398,34 @@ public static int Bisection(double x, IList values, int start = 0, SortO // Perform bisection search int xlo = start, xhi = N, xm = 0; - while (xhi - xlo > 1) + if (order == SortOrder.Ascending) { - xm = xlo + (xhi - xlo >> 1); - if (x >= values[xm] && order == SortOrder.Ascending) + while (xhi - xlo > 1) { - xlo = xm; + xm = xlo + (xhi - xlo >> 1); + if (x >= values[xm]) + { + xlo = xm; + } + else + { + xhi = xm; + } } - else + } + else + { + while (xhi - xlo > 1) { - xhi = xm; + xm = xlo + (xhi - xlo >> 1); + if (x < values[xm]) + { + xlo = xm; + } + else + { + xhi = xm; + } } } // Return XLO @@ -480,16 +498,34 @@ public static int Bisection(double x, OrderedPairedData orderedPairedData, int s // Perform bisection search int xlo = start, xhi = N, xm = 0; - while (xhi - xlo > 1) + if (orderedPairedData.OrderX == SortOrder.Ascending) { - xm = xlo + (xhi - xlo >> 1); - if (x >= orderedPairedData[xm].X && orderedPairedData.OrderX == SortOrder.Ascending) + while (xhi - xlo > 1) { - xlo = xm; + xm = xlo + (xhi - xlo >> 1); + if (x >= orderedPairedData[xm].X) + { + xlo = xm; + } + else + { + xhi = xm; + } } - else + } + else + { + while (xhi - xlo > 1) { - xhi = xm; + xm = xlo + (xhi - xlo >> 1); + if (x < orderedPairedData[xm].X) + { + xlo = xm; + } + else + { + xhi = xm; + } } } // Return XLO @@ -563,16 +599,34 @@ public static int Bisection(double x, IList ordinates, int start = 0, // Perform bisection search int xlo = start, xhi = N, xm = 0; - while (xhi - xlo > 1) + if (order == SortOrder.Ascending) { - xm = xlo + (xhi - xlo >> 1); - if (x >= ordinates[xm].X && order == SortOrder.Ascending) + while (xhi - xlo > 1) { - xlo = xm; + xm = xlo + (xhi - xlo >> 1); + if (x >= ordinates[xm].X) + { + xlo = xm; + } + else + { + xhi = xm; + } } - else + } + else + { + while (xhi - xlo > 1) { - xhi = xm; + xm = xlo + (xhi - xlo >> 1); + if (x < ordinates[xm].X) + { + xlo = xm; + } + else + { + xhi = xm; + } } } // Return XLO @@ -660,7 +714,7 @@ public static int Hunt(double x, IList values, int start = 0, SortOrder } else { - if (x >= values[xlo] && order == SortOrder.Ascending) + if ((x >= values[xlo]) == (order == SortOrder.Ascending)) { // Hunt up for (;;) @@ -668,7 +722,7 @@ public static int Hunt(double x, IList values, int start = 0, SortOrder // Not done hunting so double the increment xhi = xlo + inc; if (xhi >= N - 1) { xhi = N - 1; break; } - else if (x < values[xhi] && order == SortOrder.Ascending) break; + else if ((x < values[xhi]) == (order == SortOrder.Ascending)) break; else { xlo = xhi; @@ -684,7 +738,7 @@ public static int Hunt(double x, IList values, int start = 0, SortOrder { xlo = xlo - inc; if (xlo <= 0) { xlo = 0; break; } - else if (x >= values[xlo] && order == SortOrder.Ascending) break; + else if ((x >= values[xlo]) == (order == SortOrder.Ascending)) break; else { xhi = xlo; @@ -698,7 +752,7 @@ public static int Hunt(double x, IList values, int start = 0, SortOrder while (xhi - xlo > 1) { xm = xlo + (xhi - xlo >> 1); - if (x >= values[xm] && order == SortOrder.Ascending) + if ((x >= values[xm]) == (order == SortOrder.Ascending)) { xlo = xm; } @@ -840,7 +894,7 @@ public static int Hunt(double xValue, OrderedPairedData orderedPairedData, int s while (XHI - XLO > 1) { XM = XLO + (XHI - XLO >> 1); - if (X >= orderedPairedData[XM].X && ASCND == true) + if ((X >= orderedPairedData[XM].X) == ASCND) { XLO = XM; } @@ -988,7 +1042,7 @@ public static int Hunt(double xValue, IList ordinateData, int start = while (XHI - XLO > 1) { XM = XLO + (XHI - XLO >> 1); - if (X >= ordinateData[XM].X && ASCND == true) + if ((X >= ordinateData[XM].X) == ASCND) { XLO = XM; } diff --git a/Numerics/Data/Statistics/Histogram.cs b/Numerics/Data/Statistics/Histogram.cs index e8ec3970..52aa12b8 100644 --- a/Numerics/Data/Statistics/Histogram.cs +++ b/Numerics/Data/Statistics/Histogram.cs @@ -397,22 +397,32 @@ public void AddBin(Bin bin) public void AddData(double data) { SortBins(); - int index = GetBinIndexOf(data); if (data <= LowerBound) { - _bins.First().LowerBound = data; + if (data < LowerBound) + { + _bins.First().LowerBound = data; + LowerBound = data; + _areBinsSorted = false; + } _bins.First().Frequency += 1; - _areBinsSorted = false; } else if (data >= UpperBound) { - _bins.Last().UpperBound = data; + if (data > UpperBound) + { + _bins.Last().UpperBound = data; + UpperBound = data; + } _bins.Last().Frequency += 1; - _areBinsSorted = false; } - else if (index >= 0 && index < NumberOfBins) + else { - _bins[index].Frequency += 1; + int index = GetBinIndexOf(data); + if (index >= 0 && index < NumberOfBins) + { + _bins[index].Frequency += 1; + } } } diff --git a/Numerics/Data/Statistics/Probability.cs b/Numerics/Data/Statistics/Probability.cs index 07afaed3..6811c4ec 100644 --- a/Numerics/Data/Statistics/Probability.cs +++ b/Numerics/Data/Statistics/Probability.cs @@ -305,6 +305,7 @@ public static double JointProbabilityHPCM(IList probabilities, int[] ind // Get z-values double zMin = -9, zMax = 9; + const double minimumCdf = 1E-300; var R = new double[n, n]; Array.Copy(correlationMatrix, R, correlationMatrix.Length); int i, j, k, ir, ic; @@ -325,7 +326,7 @@ public static double JointProbabilityHPCM(IList probabilities, int[] ind z1 = R[0, 0]; pdf = Normal.StandardPDF(z1); cdf = Normal.StandardCDF(z1); - //if (cdf < 1e-300) cdf = 1e-300; + if (cdf < minimumCdf) cdf = minimumCdf; A = pdf / cdf; B = A * (z1 + A); for (k = 1; k < n; k++) @@ -352,7 +353,7 @@ public static double JointProbabilityHPCM(IList probabilities, int[] ind z1 = R[j, j - 1]; pdf = Normal.StandardPDF(z1); cdf = Normal.StandardCDF(z1); - if (cdf < 1e-300) cdf = 1e-300; + if (cdf < minimumCdf) cdf = minimumCdf; A = pdf / cdf; B = A * (z1 + A); for (k = j + 1; k < n; k++) diff --git a/Numerics/Distributions/Bivariate Copulas/AMHCopula.cs b/Numerics/Distributions/Bivariate Copulas/AMHCopula.cs index f773e28f..c236a875 100644 --- a/Numerics/Distributions/Bivariate Copulas/AMHCopula.cs +++ b/Numerics/Distributions/Bivariate Copulas/AMHCopula.cs @@ -82,6 +82,12 @@ public override double ThetaMaximum /// public override ArgumentOutOfRangeException? ValidateParameter(double parameter, bool throwException) { + if (double.IsNaN(parameter) || double.IsInfinity(parameter)) + { + var exception = new ArgumentOutOfRangeException(nameof(Theta), "The dependency parameter must be finite."); + if (throwException) throw exception; + return exception; + } if (parameter < ThetaMinimum) { if (throwException) throw new ArgumentOutOfRangeException(nameof(Theta), "The dependency parameter θ (theta) must be greater than or equal to " + ThetaMinimum.ToString() + "."); @@ -167,7 +173,7 @@ public override double[] InverseCDF(double u, double v) /// public override BivariateCopula Clone() { - return new AMHCopula(Theta, MarginalDistributionX, MarginalDistributionY); + return new AMHCopula(Theta, CloneMarginal(MarginalDistributionX), CloneMarginal(MarginalDistributionY)); } /// diff --git a/Numerics/Distributions/Bivariate Copulas/Base/ArchimedeanCopula.cs b/Numerics/Distributions/Bivariate Copulas/Base/ArchimedeanCopula.cs index b2167276..cf887913 100644 --- a/Numerics/Distributions/Bivariate Copulas/Base/ArchimedeanCopula.cs +++ b/Numerics/Distributions/Bivariate Copulas/Base/ArchimedeanCopula.cs @@ -49,6 +49,12 @@ public override void SetCopulaParameters(double[] parameters) /// public override ArgumentOutOfRangeException? ValidateParameter(double parameter, bool throwException) { + if (double.IsNaN(parameter) || double.IsInfinity(parameter)) + { + var exception = new ArgumentOutOfRangeException(nameof(Theta), "The dependency parameter must be finite."); + if (throwException) throw exception; + return exception; + } if (parameter < ThetaMinimum) { if (throwException) throw new ArgumentOutOfRangeException(nameof(Theta), "The dependency parameter θ (theta) must be greater than or equal to " + ThetaMinimum.ToString() + "."); @@ -60,7 +66,11 @@ public override void SetCopulaParameters(double[] parameters) if (throwException) throw new ArgumentOutOfRangeException(nameof(Theta), "The dependency parameter θ (theta) must be less than or equal to " + ThetaMaximum.ToString() + "."); return new ArgumentOutOfRangeException(nameof(Theta), "The dependency parameter θ (theta) must be less than or equal to " + ThetaMaximum.ToString() + "."); } - return new ArgumentOutOfRangeException(nameof(Theta),"Parameter is valid"); + // A valid parameter must return null: the Theta setter derives ParametersValid + // from 'ValidateParameter(value, false) is null', so returning a sentinel + // exception here left ParametersValid permanently false for every Archimedean + // family that did not override this method (Clayton, Gumbel, Joe). + return null; } /// diff --git a/Numerics/Distributions/Bivariate Copulas/Base/BivariateCopula.cs b/Numerics/Distributions/Bivariate Copulas/Base/BivariateCopula.cs index 17a624d5..a955fdc6 100644 --- a/Numerics/Distributions/Bivariate Copulas/Base/BivariateCopula.cs +++ b/Numerics/Distributions/Bivariate Copulas/Base/BivariateCopula.cs @@ -137,10 +137,33 @@ public double LogPDF(double u, double v) public abstract ArgumentOutOfRangeException? ValidateParameter(double parameter, bool throwException); /// - /// Create a deep copy of the copula. + /// Create a deep copy of the copula, including independently cloned marginal + /// distributions. /// + /// + /// Marginals are deep-copied through because + /// distributions memoize internal state lazily — clones sharing marginal + /// references are not safe to use concurrently (e.g., one clone per thread in + /// Monte Carlo frameworks or parallel likelihood evaluation). + /// public abstract BivariateCopula Clone(); + /// + /// Deep-copies an attached marginal distribution for . + /// + /// The marginal distribution, or null when unattached. + /// An independent copy of the marginal; null when unattached. A marginal + /// implementation outside (which carries the + /// clone support) is returned by reference, preserving the previous behavior. + protected static IUnivariateDistribution? CloneMarginal(IUnivariateDistribution? marginal) + { + if (marginal is UnivariateDistributionBase cloneable) + { + return cloneable.Clone(); + } + return marginal; + } + /// /// Returns the OR joint exceedance probability. When either of the variables exceeds a particular threshold value /// diff --git a/Numerics/Distributions/Bivariate Copulas/ClaytonCopula.cs b/Numerics/Distributions/Bivariate Copulas/ClaytonCopula.cs index 0ce2edbc..ade724e6 100644 --- a/Numerics/Distributions/Bivariate Copulas/ClaytonCopula.cs +++ b/Numerics/Distributions/Bivariate Copulas/ClaytonCopula.cs @@ -151,7 +151,7 @@ public override double LowerTailDependence /// public override BivariateCopula Clone() { - return new ClaytonCopula(Theta, MarginalDistributionX, MarginalDistributionY); + return new ClaytonCopula(Theta, CloneMarginal(MarginalDistributionX), CloneMarginal(MarginalDistributionY)); } /// diff --git a/Numerics/Distributions/Bivariate Copulas/FrankCopula.cs b/Numerics/Distributions/Bivariate Copulas/FrankCopula.cs index 1a976edb..15868b7f 100644 --- a/Numerics/Distributions/Bivariate Copulas/FrankCopula.cs +++ b/Numerics/Distributions/Bivariate Copulas/FrankCopula.cs @@ -81,6 +81,12 @@ public override double ThetaMaximum /// public override ArgumentOutOfRangeException? ValidateParameter(double parameter, bool throwException) { + if (double.IsNaN(parameter) || double.IsInfinity(parameter)) + { + var exception = new ArgumentOutOfRangeException(nameof(Theta), "The dependency parameter must be finite."); + if (throwException) throw exception; + return exception; + } if (parameter < ThetaMinimum) { if (throwException) throw new ArgumentOutOfRangeException(nameof(Theta), "The dependency parameter θ (theta) must be greater than or equal to " + ThetaMinimum.ToString() + "."); @@ -169,7 +175,7 @@ public override double[] InverseCDF(double u, double v) /// public override BivariateCopula Clone() { - return new FrankCopula(Theta, MarginalDistributionX, MarginalDistributionY); + return new FrankCopula(Theta, CloneMarginal(MarginalDistributionX), CloneMarginal(MarginalDistributionY)); } /// diff --git a/Numerics/Distributions/Bivariate Copulas/GumbelCopula.cs b/Numerics/Distributions/Bivariate Copulas/GumbelCopula.cs index f6819fbc..d2bd32d5 100644 --- a/Numerics/Distributions/Bivariate Copulas/GumbelCopula.cs +++ b/Numerics/Distributions/Bivariate Copulas/GumbelCopula.cs @@ -149,7 +149,7 @@ public override double UpperTailDependence /// public override BivariateCopula Clone() { - return new GumbelCopula(Theta, MarginalDistributionX, MarginalDistributionY); + return new GumbelCopula(Theta, CloneMarginal(MarginalDistributionX), CloneMarginal(MarginalDistributionY)); } /// diff --git a/Numerics/Distributions/Bivariate Copulas/JoeCopula.cs b/Numerics/Distributions/Bivariate Copulas/JoeCopula.cs index 1d2c66c9..32fe731e 100644 --- a/Numerics/Distributions/Bivariate Copulas/JoeCopula.cs +++ b/Numerics/Distributions/Bivariate Copulas/JoeCopula.cs @@ -151,7 +151,7 @@ public override double UpperTailDependence /// public override BivariateCopula Clone() { - return new JoeCopula(Theta, MarginalDistributionX, MarginalDistributionY); + return new JoeCopula(Theta, CloneMarginal(MarginalDistributionX), CloneMarginal(MarginalDistributionY)); } /// diff --git a/Numerics/Distributions/Bivariate Copulas/NormalCopula.cs b/Numerics/Distributions/Bivariate Copulas/NormalCopula.cs index 4ca6bb3e..9adf1991 100644 --- a/Numerics/Distributions/Bivariate Copulas/NormalCopula.cs +++ b/Numerics/Distributions/Bivariate Copulas/NormalCopula.cs @@ -97,6 +97,12 @@ public override double ThetaMaximum /// public override ArgumentOutOfRangeException? ValidateParameter(double parameter, bool throwException) { + if (double.IsNaN(parameter) || double.IsInfinity(parameter)) + { + var exception = new ArgumentOutOfRangeException(nameof(Theta), "The correlation parameter must be finite."); + if (throwException) throw exception; + return exception; + } if (parameter < ThetaMinimum) { if (throwException) throw new ArgumentOutOfRangeException(nameof(Theta), "The correlation parameter ρ (rho) must be greater than " + ThetaMinimum.ToString() + "."); @@ -178,7 +184,7 @@ public override double[] InverseCDF(double u, double v) /// public override BivariateCopula Clone() { - return new NormalCopula(Theta, MarginalDistributionX, MarginalDistributionY); + return new NormalCopula(Theta, CloneMarginal(MarginalDistributionX), CloneMarginal(MarginalDistributionY)); } } diff --git a/Numerics/Distributions/Bivariate Copulas/StudentTCopula.cs b/Numerics/Distributions/Bivariate Copulas/StudentTCopula.cs index f2fdc36b..4235229d 100644 --- a/Numerics/Distributions/Bivariate Copulas/StudentTCopula.cs +++ b/Numerics/Distributions/Bivariate Copulas/StudentTCopula.cs @@ -189,6 +189,12 @@ public override double ThetaMaximum /// When true, throws on the first invalid parameter; when false, returns the exception instead. public ArgumentOutOfRangeException? ValidateParameters(double rho, double degreesOfFreedom, bool throwException) { + if (double.IsNaN(rho) || double.IsInfinity(rho)) + { + var exception = new ArgumentOutOfRangeException(nameof(Theta), "The correlation parameter must be finite."); + if (throwException) throw exception; + return exception; + } if (rho < ThetaMinimum) { if (throwException) throw new ArgumentOutOfRangeException(nameof(Theta), "The correlation parameter ρ (rho) must be greater than " + ThetaMinimum.ToString() + "."); @@ -409,7 +415,7 @@ public override double UpperTailDependence /// public override BivariateCopula Clone() { - return new StudentTCopula(Theta, _nu, MarginalDistributionX, MarginalDistributionY); + return new StudentTCopula(Theta, _nu, CloneMarginal(MarginalDistributionX), CloneMarginal(MarginalDistributionY)); } } diff --git a/Numerics/Distributions/Multivariate/BivariateEmpirical.cs b/Numerics/Distributions/Multivariate/BivariateEmpirical.cs index c9b92e57..018628d4 100644 --- a/Numerics/Distributions/Multivariate/BivariateEmpirical.cs +++ b/Numerics/Distributions/Multivariate/BivariateEmpirical.cs @@ -143,6 +143,7 @@ public void SetParameters(IList x1Values, IList x2Values, double X1Values = x1Values.ToArray(); X2Values = x2Values.ToArray(); ProbabilityValues = pValues; + bilinear = null; } @@ -170,6 +171,14 @@ public void SetParameters(IList x1Values, IList x2Values, double } // Check if X1 and X2 are in ascending order + for (int i = 0; i < x1Values.Count; i++) + { + if (double.IsNaN(x1Values[i]) || double.IsInfinity(x1Values[i])) + { + if (throwException) throw new ArgumentOutOfRangeException("Primary Values", "Primary values must be finite."); + return new ArgumentOutOfRangeException("Primary Values", "Primary values must be finite."); + } + } for (int i = 1; i < x1Values.Count; i++) { if (x1Values[i] <= x1Values[i - 1]) @@ -179,6 +188,14 @@ public void SetParameters(IList x1Values, IList x2Values, double return new ArgumentOutOfRangeException("Primary Values", "Primary values must be in ascending order."); } } + for (int i = 0; i < x2Values.Count; i++) + { + if (double.IsNaN(x2Values[i]) || double.IsInfinity(x2Values[i])) + { + if (throwException) throw new ArgumentOutOfRangeException("Secondary Values", "Secondary values must be finite."); + return new ArgumentOutOfRangeException("Secondary Values", "Secondary values must be finite."); + } + } for (int i = 1; i < x2Values.Count; i++) { if (x2Values[i] <= x2Values[i - 1]) @@ -207,7 +224,7 @@ public void SetParameters(IList x1Values, IList x2Values, double { for (int j = 0; j < pValues.GetLength(1); j++) // Columns { - if (pValues[i, j] < 0.0d || pValues[i, j] > 1.0d) + if (double.IsNaN(pValues[i, j]) || double.IsInfinity(pValues[i, j]) || pValues[i, j] < 0.0d || pValues[i, j] > 1.0d) { if (throwException) throw new ArgumentOutOfRangeException(nameof(ProbabilityValues), "Probability values must be equal to or between 0 and 1."); return new ArgumentOutOfRangeException(nameof(ProbabilityValues), "Probability values must be equal to or between 0 and 1."); diff --git a/Numerics/Distributions/Multivariate/MultivariateNormal.cs b/Numerics/Distributions/Multivariate/MultivariateNormal.cs index 90cedc82..b0d97929 100644 --- a/Numerics/Distributions/Multivariate/MultivariateNormal.cs +++ b/Numerics/Distributions/Multivariate/MultivariateNormal.cs @@ -300,7 +300,15 @@ private void CreateCorrelationMatrix() { var ex = new ArgumentOutOfRangeException(nameof(covariance), "Covariance matrix must not be null."); if (throwException) throw ex; else return ex; - } + } + for (int i = 0; i < mean.Length; i++) + { + if (double.IsNaN(mean[i]) || double.IsInfinity(mean[i])) + { + var ex = new ArgumentOutOfRangeException(nameof(mean), "Mean values must be finite."); + if (throwException) throw ex; else return ex; + } + } var m = new Matrix(covariance); if (!m.IsSquare) { @@ -312,6 +320,17 @@ private void CreateCorrelationMatrix() var ex = new ArgumentOutOfRangeException(nameof(Covariance), "Mean length must match covariance dimension."); if (throwException) throw ex; else return ex; } + for (int i = 0; i < m.NumberOfRows; i++) + { + for (int j = 0; j < m.NumberOfColumns; j++) + { + if (double.IsNaN(m[i, j]) || double.IsInfinity(m[i, j])) + { + var ex = new ArgumentOutOfRangeException(nameof(covariance), "Covariance values must be finite."); + if (throwException) throw ex; else return ex; + } + } + } var chol = new CholeskyDecomposition(m); if (!chol.IsPositiveDefinite) @@ -322,12 +341,199 @@ private void CreateCorrelationMatrix() return null; } + /// + /// Attempts to set the distribution parameters without throwing: the covariance + /// is factorized eagerly, and when it is not positive-definite the density is + /// marked invalid — returns negative infinity and + /// returns zero until a valid covariance is set. + /// + /// The mean vector μ (mu) for the distribution. + /// The covariance matrix Σ (sigma) for the distribution. + /// True when the covariance is positive-definite and the density is usable. + /// + /// The non-throwing mutable path for samplers and likelihood loops that swap + /// covariances per evaluation (e.g., MCMC proposals of correlation parameters): + /// an infeasible proposal scores −∞ and is rejected instead of raising an + /// exception. keeps its throwing contract. + /// + public bool TrySetParameters(double[] mean, double[,] covariance) + { + // The Cholesky constructor itself throws on strongly indefinite matrices + // (negative pivots) while merely flagging weakly non-positive-definite + // ones, so the non-throwing contract absorbs both failure modes. + try + { + if (ValidateParameters(mean, covariance, false) is null) + { + SetParameters(mean, covariance); + _densityValid = true; + return true; + } + } + catch (Exception ex) + { + System.Diagnostics.Debug.WriteLine( + $"TrySetParameters rejected a covariance: {ex.Message}"); + } + _densityValid = false; + return false; + } + + /// + /// Attempts to swap the covariance matrix in place (the mean is kept) without + /// throwing; see for the invalid-density + /// contract. + /// + /// The covariance matrix Σ (sigma) for the distribution. + /// True when the covariance is positive-definite and the density is usable. + public bool TrySetCovariance(double[,] covariance) + { + return TrySetParameters(_mean, covariance); + } + + /// + /// Whether the current covariance factorization supports density evaluation. + /// False after a failed / + /// until a valid covariance is set. + /// + public bool IsDensityValid => _densityValid; + + /// Tracks the non-throwing mutable-covariance state. + private bool _densityValid = true; + + /// + /// Returns the marginal distribution over a subset of the dimensions — the + /// sub-mean and sub-covariance at the given indices (the marginal of a joint + /// Gaussian is the corresponding sub-Gaussian). + /// + /// The zero-based dimension indices to keep, in the output order. + /// The marginal distribution. + /// Thrown when the indices are null, empty, repeated, or out of range. + public MultivariateNormal Marginal(params int[] indices) + { + ValidateIndices(indices, Dimension); + var mean = new double[indices.Length]; + var covariance = new double[indices.Length, indices.Length]; + for (int i = 0; i < indices.Length; i++) + { + mean[i] = _mean[indices[i]]; + for (int j = 0; j < indices.Length; j++) + { + covariance[i, j] = _covariance[indices[i], indices[j]]; + } + } + return new MultivariateNormal(mean, covariance); + } + + /// + /// Returns the conditional distribution of the remaining dimensions given + /// observed values at a subset of the dimensions: + /// μ_c = μ₁ + Σ₁₂·Σ₂₂⁻¹·(x₂ − μ₂) and Σ_c = Σ₁₁ − Σ₁₂·Σ₂₂⁻¹·Σ₂₁ (the Schur + /// complement), with Σ₂₂ solved through its Cholesky factorization. + /// + /// The zero-based dimension indices of the observed subset. + /// The observed values, parallel to the indices. + /// The conditional distribution over the remaining dimensions, in ascending dimension order. + /// Thrown when the indices are invalid, the values length differs, or every dimension is observed. + /// + /// Reference: Eaton, M.L. (1983). Multivariate Statistics: A Vector Space + /// Approach. Wiley. (The Gaussian conditioning identities.) + /// + public MultivariateNormal Conditional(int[] observedIndices, double[] observedValues) + { + ValidateIndices(observedIndices, Dimension); + if (observedValues == null || observedValues.Length != observedIndices.Length) + throw new ArgumentOutOfRangeException(nameof(observedValues), "The observed values must parallel the observed indices."); + if (observedIndices.Length >= Dimension) + throw new ArgumentOutOfRangeException(nameof(observedIndices), "At least one dimension must remain unobserved."); + + // The remaining (free) dimensions in ascending order. + var isObserved = new bool[Dimension]; + foreach (int index in observedIndices) + isObserved[index] = true; + var free = new int[Dimension - observedIndices.Length]; + int f = 0; + for (int i = 0; i < Dimension; i++) + { + if (!isObserved[i]) free[f++] = i; + } + + // Partition: Σ22 (observed), Σ12 (free × observed), and the innovation + // z = x₂ − μ₂. + int nObserved = observedIndices.Length, nFree = free.Length; + var sigma22 = new Matrix(nObserved, nObserved); + for (int i = 0; i < nObserved; i++) + { + for (int j = 0; j < nObserved; j++) + sigma22[i, j] = _covariance[observedIndices[i], observedIndices[j]]; + } + var cholesky22 = new CholeskyDecomposition(sigma22); + if (!cholesky22.IsPositiveDefinite) + throw new ArgumentOutOfRangeException(nameof(observedIndices), "The observed sub-covariance is not positive-definite."); + + var z = new double[nObserved]; + for (int i = 0; i < nObserved; i++) + z[i] = observedValues[i] - _mean[observedIndices[i]]; + var solvedZ = cholesky22.Solve(new Vector(z)); + + // Solve Σ22⁻¹·Σ21 column-by-column (Σ21 columns are the free rows of Σ12ᵀ). + var solvedColumns = new Vector[nFree]; + for (int i = 0; i < nFree; i++) + { + var column = new double[nObserved]; + for (int j = 0; j < nObserved; j++) + column[j] = _covariance[observedIndices[j], free[i]]; + solvedColumns[i] = cholesky22.Solve(new Vector(column)); + } + + var conditionalMean = new double[nFree]; + var conditionalCovariance = new double[nFree, nFree]; + for (int i = 0; i < nFree; i++) + { + double meanAdjustment = 0d; + for (int j = 0; j < nObserved; j++) + meanAdjustment += _covariance[free[i], observedIndices[j]] * solvedZ[j]; + conditionalMean[i] = _mean[free[i]] + meanAdjustment; + for (int k = 0; k < nFree; k++) + { + double schur = 0d; + for (int j = 0; j < nObserved; j++) + schur += _covariance[free[i], observedIndices[j]] * solvedColumns[k][j]; + conditionalCovariance[i, k] = _covariance[free[i], free[k]] - schur; + } + } + return new MultivariateNormal(conditionalMean, conditionalCovariance); + } + + /// + /// Validates a dimension-index subset: non-null, non-empty, in range, and + /// without repeats. + /// + /// The indices to validate. + /// The distribution dimension. + /// Thrown when the subset is invalid. + private static void ValidateIndices(int[] indices, int dimension) + { + if (indices == null || indices.Length == 0) + throw new ArgumentOutOfRangeException(nameof(indices), "At least one dimension index is required."); + var seen = new bool[dimension]; + foreach (int index in indices) + { + if (index < 0 || index >= dimension) + throw new ArgumentOutOfRangeException(nameof(indices), "A dimension index is out of range."); + if (seen[index]) + throw new ArgumentOutOfRangeException(nameof(indices), "Dimension indices must not repeat."); + seen[index] = true; + } + } + /// /// The Probability Density Function (PDF) of the distribution evaluated at a point X. /// /// A point in the distribution space. public override double PDF(double[] x) { + if (!_densityValid) return 0d; return Math.Exp(-0.5d * Mahalanobis(x) + _lnconstant); } @@ -337,6 +543,7 @@ public override double PDF(double[] x) /// The vector of x values. public override double LogPDF(double[] x) { + if (!_densityValid) return double.NegativeInfinity; double f = -0.5d * Mahalanobis(x) + _lnconstant; if (double.IsNaN(f) || double.IsInfinity(f)) return double.NegativeInfinity; return f; @@ -1188,7 +1395,7 @@ public static double BivariateCDF(double z1, double z2, double r) /// array of lower integration limits. /// array of upper integration limits. /// array of integration limits flags: - /// if INFIN(I) > 0, Ith limits are (-infinity, infinity); + /// if INFIN(I) < 0, Ith limits are (-infinity, infinity); /// if INFIN(I) = 0, Ith limits are(-infinity, UPPER(I)]; /// if INFIN(I) = 1, Ith limits are[LOWER(I), infinity); /// if INFIN(I) = 2, Ith limits are[LOWER(I), UPPER(I)]. @@ -1214,7 +1421,8 @@ public void MVNDST(int N, double[] LOWER, double[] UPPER, int[] INFIN, double[] } else { - INFORM = (int)MVNDNT(N, CORREL, LOWER, UPPER, INFIN, ref INFIS, ref D, ref E, Y); + INFORM = 0; + MVNDNT(N, CORREL, LOWER, UPPER, INFIN, ref INFIS, ref D, ref E, Y); if (N - INFIS == 0) { VALUE = 1; @@ -1311,10 +1519,25 @@ private double MVNDFN(int N, double[] W) return result; } - private double MVNDNT(int N, double[] CORREL, double[] LOWER, double[] UPPER, int[] INFIN, ref int INFIS, ref double D, ref double E, double[] Y) + /// + /// Initializes the transformed limits and covariance factor used by the multivariate-normal integrand. + /// + /// The number of integration dimensions. + /// The packed strict-lower-triangle correlation coefficients. + /// The lower integration limits. + /// The upper integration limits. + /// The integration-bound flags. + /// Receives the number of dimensions having no finite bound, including any analytically collapsed dimension. + /// Receives the lower transformed probability for an analytically evaluated one- or two-dimensional problem. + /// Receives the upper transformed probability for an analytically evaluated one- or two-dimensional problem. + /// The conditioning work vector used while sorting the covariance matrix. + /// + /// This method is the C# translation of the MVNDNT initialization entry point in Alan Genz's + /// Fortran MVNDFN routine. The original function returns zero solely to initialize the caller's + /// INFORM status; performs that initialization explicitly in C#. + /// + private void MVNDNT(int N, double[] CORREL, double[] LOWER, double[] UPPER, int[] INFIN, ref int INFIS, ref double D, ref double E, double[] Y) { - double result = 0; - COVSRT(N, LOWER, UPPER, CORREL, INFIN, Y, ref INFIS, MVNDFN_A, MVNDFN_B, MVNDFN_COV, MVNDFN_INFI); if (N - INFIS == 1) @@ -1348,8 +1571,6 @@ private double MVNDNT(int N, double[] CORREL, double[] LOWER, double[] UPPER, in } INFIS++; } - - return result; } private void MVNLMS(double A, double B, int INFIN, ref double LOWER, ref double UPPER) @@ -1367,17 +1588,21 @@ private void MVNLMS(double A, double B, int INFIN, ref double LOWER, ref double /// /// Subroutine to sort integration limits and determine Cholesky factor. /// - /// - /// - /// - /// - /// - /// - /// - /// - /// - /// - /// + /// The number of integration dimensions. + /// The lower integration limits, reordered in place. + /// The upper integration limits, reordered in place. + /// The packed strict-lower-triangle correlation coefficients. + /// The bound flags: 0 for a lower-infinite interval, 1 for an upper-infinite interval, and 2 for finite bounds. + /// The conditioning work vector. + /// Receives the number of dimensions with no finite bound. + /// Receives the sorted and standardized lower limits. + /// Receives the sorted and standardized upper limits. + /// Receives the packed lower-triangular covariance factor. + /// Receives the sorted bound flags. + /// + /// This is a zero-based translation of Alan Genz's MVNDST covariance sorting routine. + /// The degenerate branch preserves the packed-array permutation used for singular covariance matrices. + /// private void COVSRT(int N, double[] LOWER, double[] UPPER, double[] CORREL, int[] INFIN, double[] Y, ref int INFIS, double[] A, double[] B, double[] COV, int[] INFI) { double SQTWPI = 2.506628274631001; @@ -1535,7 +1760,7 @@ private void COVSRT(int N, double[] LOWER, double[] UPPER, double[] CORREL, int[ // If the covariance matrix diagonal entry is zero, permute limits and/ or rows, if necessary. - for (int j = i - 1; j < 0; j--) + for (int j = i - 1; j >= 0; j--) { if (Math.Abs(COV[II + j]) > EPS) { @@ -1546,27 +1771,27 @@ private void COVSRT(int N, double[] LOWER, double[] UPPER, double[] CORREL, int[ Tools.Swap(ref A[i], ref B[i]); if (INFI[i] != 2) INFI[i] = 1 - INFI[i]; } - for (int l = 0; l < j; l++) + for (int l = 0; l <= j; l++) { COV[II + l] = COV[II + l] / COV[II + j]; } - for (int l = j + 1; l < i - 1; l++) + for (int l = j + 1; l <= i - 1; l++) { - if (COV[(l - 1) * l / 2 + j + 1] > 0) + if (COV[l * (l + 1) / 2 + j + 1] > 0) { IJ = II; - for (int k = i - 1; k < l; k--) + for (int k = i - 1; k >= l; k--) { - for (int m = 0; m < k; m++) + for (int m = 0; m <= k; m++) { - Tools.Swap(ref COV[IJ - k + m], ref COV[IJ + m]); + Tools.Swap(ref COV[IJ - k + m - 1], ref COV[IJ + m]); } Tools.Swap(ref A[k], ref A[k + 1]); Tools.Swap(ref B[k], ref B[k + 1]); M = INFI[k]; INFI[k] = INFI[k + 1]; INFI[k + 1] = M; - IJ = IJ - k; + IJ = IJ - k - 1; } goto twenty; } diff --git a/Numerics/Distributions/Multivariate/MultivariateStudentT.cs b/Numerics/Distributions/Multivariate/MultivariateStudentT.cs index 5282f53a..df61a289 100644 --- a/Numerics/Distributions/Multivariate/MultivariateStudentT.cs +++ b/Numerics/Distributions/Multivariate/MultivariateStudentT.cs @@ -314,7 +314,7 @@ public void SetParameters(double degreesOfFreedom, double[] location, double[,] /// public ArgumentOutOfRangeException? ValidateParameters(double degreesOfFreedom, double[] location, double[,] scaleMatrix, bool throwException) { - if (degreesOfFreedom <= 0) + if (double.IsNaN(degreesOfFreedom) || double.IsInfinity(degreesOfFreedom) || degreesOfFreedom <= 0) { var ex = new ArgumentOutOfRangeException(nameof(degreesOfFreedom), "Degrees of freedom must be greater than zero."); if (throwException) throw ex; else return ex; @@ -329,6 +329,14 @@ public void SetParameters(double degreesOfFreedom, double[] location, double[,] var ex = new ArgumentOutOfRangeException(nameof(scaleMatrix), "Scale matrix must not be null."); if (throwException) throw ex; else return ex; } + for (int i = 0; i < location.Length; i++) + { + if (double.IsNaN(location[i]) || double.IsInfinity(location[i])) + { + var ex = new ArgumentOutOfRangeException(nameof(location), "Location values must be finite."); + if (throwException) throw ex; else return ex; + } + } var m = new Matrix(scaleMatrix); if (!m.IsSquare) { @@ -340,6 +348,17 @@ public void SetParameters(double degreesOfFreedom, double[] location, double[,] var ex = new ArgumentOutOfRangeException(nameof(scaleMatrix), "Location vector length must match scale matrix dimension."); if (throwException) throw ex; else return ex; } + for (int i = 0; i < m.NumberOfRows; i++) + { + for (int j = 0; j < m.NumberOfColumns; j++) + { + if (double.IsNaN(m[i, j]) || double.IsInfinity(m[i, j])) + { + var ex = new ArgumentOutOfRangeException(nameof(scaleMatrix), "Scale-matrix values must be finite."); + if (throwException) throw ex; else return ex; + } + } + } try { var chol = new CholeskyDecomposition(m); diff --git a/Numerics/Distributions/Univariate/Base/UnivariateDistributionFactory.cs b/Numerics/Distributions/Univariate/Base/UnivariateDistributionFactory.cs index 58cad0b5..b9179f88 100644 --- a/Numerics/Distributions/Univariate/Base/UnivariateDistributionFactory.cs +++ b/Numerics/Distributions/Univariate/Base/UnivariateDistributionFactory.cs @@ -23,155 +23,105 @@ public sealed class UnivariateDistributionFactory /// /// A univariate distribution. /// + /// + /// The requested distribution type requires additional component distributions or a user implementation. + /// + /// + /// is not a defined value. + /// public static UnivariateDistributionBase CreateDistribution(UnivariateDistributionType distributionType) { - if (distributionType == UnivariateDistributionType.Bernoulli) - { - return new Bernoulli(); - } - else if (distributionType == UnivariateDistributionType.Beta) - { - return new BetaDistribution(); - } - else if (distributionType == UnivariateDistributionType.Binomial) - { - return new Binomial(); - } - else if (distributionType == UnivariateDistributionType.Cauchy) - { - return new Cauchy(); - } - else if (distributionType == UnivariateDistributionType.ChiSquared) - { - return new ChiSquared(); - } - else if (distributionType == UnivariateDistributionType.Exponential) - { - return new Exponential(); - } - else if (distributionType == UnivariateDistributionType.GammaDistribution) - { - return new GammaDistribution(); - } - else if (distributionType == UnivariateDistributionType.GeneralizedBeta) - { - return new GeneralizedBeta(); - } - else if (distributionType == UnivariateDistributionType.GeneralizedExtremeValue) - { - return new GeneralizedExtremeValue(); - } - else if (distributionType == UnivariateDistributionType.GeneralizedLogistic) - { - return new GeneralizedLogistic(); - } - else if (distributionType == UnivariateDistributionType.GeneralizedNormal) - { - return new GeneralizedNormal(); - } - else if (distributionType == UnivariateDistributionType.GeneralizedPareto) - { - return new GeneralizedPareto(); - } - else if (distributionType == UnivariateDistributionType.Geometric) - { - return new Geometric(); - } - else if (distributionType == UnivariateDistributionType.Gumbel) - { - return new Gumbel(); - } - else if (distributionType == UnivariateDistributionType.InverseChiSquared) - { - return new InverseChiSquared(); - } - else if (distributionType == UnivariateDistributionType.InverseGamma) - { - return new InverseGamma(); - } - else if (distributionType == UnivariateDistributionType.KappaFour) - { - return new KappaFour(); - } - else if (distributionType == UnivariateDistributionType.LnNormal) - { - return new LnNormal(); - } - else if (distributionType == UnivariateDistributionType.Logistic) - { - return new Logistic(); - } - else if (distributionType == UnivariateDistributionType.LogNormal) - { - return new LogNormal(); - } - else if (distributionType == UnivariateDistributionType.LogPearsonTypeIII) - { - return new LogPearsonTypeIII(); - } - else if (distributionType == UnivariateDistributionType.NoncentralT) - { - return new NoncentralT(); - } - else if (distributionType == UnivariateDistributionType.Normal) - { - return new Normal(); - } - else if (distributionType == UnivariateDistributionType.Pareto) - { - return new Pareto(); - } - else if (distributionType == UnivariateDistributionType.PearsonTypeIII) - { - return new PearsonTypeIII(); - } - else if (distributionType == UnivariateDistributionType.Pert) - { - return new Pert(); - } - else if (distributionType == UnivariateDistributionType.PertPercentile) - { - return new PertPercentile(); - } - else if (distributionType == UnivariateDistributionType.PertPercentileZ) - { - return new PertPercentileZ(); - } - else if (distributionType == UnivariateDistributionType.Poisson) - { - return new Poisson(); - } - else if (distributionType == UnivariateDistributionType.Rayleigh) - { - return new Rayleigh(); - } - else if (distributionType == UnivariateDistributionType.StudentT) - { - return new StudentT(); - } - else if (distributionType == UnivariateDistributionType.Triangular) - { - return new Triangular(); + switch (distributionType) + { + case UnivariateDistributionType.ChiSquared: + return new ChiSquared(); + case UnivariateDistributionType.Bernoulli: + return new Bernoulli(); + case UnivariateDistributionType.Beta: + return new BetaDistribution(); + case UnivariateDistributionType.Binomial: + return new Binomial(); + case UnivariateDistributionType.Cauchy: + return new Cauchy(); + case UnivariateDistributionType.Deterministic: + return new Deterministic(); + case UnivariateDistributionType.Empirical: + return new EmpiricalDistribution(); + case UnivariateDistributionType.Exponential: + return new Exponential(); + case UnivariateDistributionType.GammaDistribution: + return new GammaDistribution(); + case UnivariateDistributionType.GeneralizedBeta: + return new GeneralizedBeta(); + case UnivariateDistributionType.GeneralizedExtremeValue: + return new GeneralizedExtremeValue(); + case UnivariateDistributionType.GeneralizedLogistic: + return new GeneralizedLogistic(); + case UnivariateDistributionType.GeneralizedNormal: + return new GeneralizedNormal(); + case UnivariateDistributionType.GeneralizedPareto: + return new GeneralizedPareto(); + case UnivariateDistributionType.Geometric: + return new Geometric(); + case UnivariateDistributionType.Gumbel: + return new Gumbel(); + case UnivariateDistributionType.InverseChiSquared: + return new InverseChiSquared(); + case UnivariateDistributionType.InverseGamma: + return new InverseGamma(); + case UnivariateDistributionType.KappaFour: + return new KappaFour(); + case UnivariateDistributionType.KernelDensity: + return new KernelDensity(); + case UnivariateDistributionType.LnNormal: + return new LnNormal(); + case UnivariateDistributionType.Logistic: + return new Logistic(); + case UnivariateDistributionType.LogNormal: + return new LogNormal(); + case UnivariateDistributionType.LogPearsonTypeIII: + return new LogPearsonTypeIII(); + case UnivariateDistributionType.NoncentralT: + return new NoncentralT(); + case UnivariateDistributionType.Normal: + return new Normal(); + case UnivariateDistributionType.Pareto: + return new Pareto(); + case UnivariateDistributionType.PearsonTypeIII: + return new PearsonTypeIII(); + case UnivariateDistributionType.Pert: + return new Pert(); + case UnivariateDistributionType.PertPercentile: + return new PertPercentile(); + case UnivariateDistributionType.PertPercentileZ: + return new PertPercentileZ(); + case UnivariateDistributionType.Poisson: + return new Poisson(); + case UnivariateDistributionType.Rayleigh: + return new Rayleigh(); + case UnivariateDistributionType.StudentT: + return new StudentT(); + case UnivariateDistributionType.Triangular: + return new Triangular(); + case UnivariateDistributionType.TruncatedNormal: + return new TruncatedNormal(); + case UnivariateDistributionType.Uniform: + return new Uniform(); + case UnivariateDistributionType.UniformDiscrete: + return new UniformDiscrete(); + case UnivariateDistributionType.VonMises: + return new VonMises(); + case UnivariateDistributionType.Weibull: + return new Weibull(); + case UnivariateDistributionType.CompetingRisks: + case UnivariateDistributionType.Mixture: + case UnivariateDistributionType.UserDefined: + throw new NotSupportedException( + "Distribution type " + distributionType + " requires external components and cannot be created without parameters."); + default: + throw new ArgumentOutOfRangeException( + nameof(distributionType), distributionType, "The distribution type is not defined."); } - else if (distributionType == UnivariateDistributionType.TruncatedNormal) - { - return new TruncatedNormal(); - } - else if (distributionType == UnivariateDistributionType.Uniform) - { - return new Uniform(); - } - else if (distributionType == UnivariateDistributionType.UniformDiscrete) - { - return new UniformDiscrete(); - } - else if (distributionType == UnivariateDistributionType.Weibull) - { - return new Weibull(); - } - - // Default to Deterministic for unrecognized types - return new Deterministic(); } /// @@ -181,6 +131,12 @@ public static UnivariateDistributionBase CreateDistribution(UnivariateDistributi /// /// A univariate distribution. /// + /// + /// The serialized distribution type requires a user-provided implementation. + /// + /// + /// The serialized distribution type is not a defined value. + /// public static UnivariateDistributionBase CreateDistribution(XElement xElement) { UnivariateDistributionType type = UnivariateDistributionType.Deterministic; diff --git a/Numerics/Distributions/Univariate/BetaDistribution.cs b/Numerics/Distributions/Univariate/BetaDistribution.cs index c09f3dc2..08922353 100644 --- a/Numerics/Distributions/Univariate/BetaDistribution.cs +++ b/Numerics/Distributions/Univariate/BetaDistribution.cs @@ -146,9 +146,11 @@ public override double Mode { if (Alpha == 1d && Beta == 1d) return 0.5d; - if (Alpha <= 1d && Beta > 1d) + if (Alpha < 1d && Beta < 1d) + return 0.5d; + if (Alpha <= 1d && Beta >= 1d) return 0.0d; - if (Alpha > 1d && Beta <= 1d) + if (Alpha >= 1d && Beta <= 1d) return 1.0d; return (Alpha - 1d) / (Alpha + Beta - 2d); } diff --git a/Numerics/Distributions/Univariate/Binomial.cs b/Numerics/Distributions/Univariate/Binomial.cs index c5073359..fbc295f9 100644 --- a/Numerics/Distributions/Univariate/Binomial.cs +++ b/Numerics/Distributions/Univariate/Binomial.cs @@ -206,8 +206,9 @@ public override double[] MaximumOfParameters /// public override void SetParameters(IList parameters) { - ProbabilityOfSuccess = parameters[0]; - NumberOfTrials = (int)parameters[1]; + _probabilityOfSuccess = parameters[0]; + _numberOfTrials = (int)parameters[1]; + _parametersValid = ValidateParameters(parameters, false) is null; } /// diff --git a/Numerics/Distributions/Univariate/ChiSquared.cs b/Numerics/Distributions/Univariate/ChiSquared.cs index f2b13531..b5a74a76 100644 --- a/Numerics/Distributions/Univariate/ChiSquared.cs +++ b/Numerics/Distributions/Univariate/ChiSquared.cs @@ -189,7 +189,8 @@ public override double[] MaximumOfParameters /// The degrees of freedom ν (nu). Range: ν > 0. public void SetParameters(double v) { - DegreesOfFreedom = (int)v; + _degreesOfFreedom = (int)v; + _parametersValid = !double.IsNaN(v) && !double.IsInfinity(v) && ValidateParameters(_degreesOfFreedom, false) is null; } /// @@ -218,6 +219,12 @@ public override void SetParameters(IList parameters) /// public override ArgumentOutOfRangeException? ValidateParameters(IList parameters, bool throwException) { + if (double.IsNaN(parameters[0]) || double.IsInfinity(parameters[0])) + { + var exception = new ArgumentOutOfRangeException(nameof(DegreesOfFreedom), "The degrees of freedom must be finite."); + if (throwException) throw exception; + return exception; + } return ValidateParameters((int)parameters[0], throwException); } diff --git a/Numerics/Distributions/Univariate/CompetingRisks.cs b/Numerics/Distributions/Univariate/CompetingRisks.cs index a963b188..9416feff 100644 --- a/Numerics/Distributions/Univariate/CompetingRisks.cs +++ b/Numerics/Distributions/Univariate/CompetingRisks.cs @@ -56,6 +56,7 @@ public CompetingRisks(IUnivariateDistribution[] distributions) private bool _empiricalCDFCreated = false; private double[,] _correlationMatrix = null!; private bool _mvnCreated = false; + private Probability.DependencyType _dependency = Probability.DependencyType.Independent; private MultivariateNormal _mvn = null!; // Soft finite floor used in tail arithmetic before returning the final log-density. @@ -85,7 +86,18 @@ public CompetingRisks(IUnivariateDistribution[] distributions) /// /// The dependency between random variables. /// - public Probability.DependencyType Dependency { get; set; } = Probability.DependencyType.Independent; + public Probability.DependencyType Dependency + { + get { return _dependency; } + set + { + if (_dependency != value) + { + _dependency = value; + _mvnCreated = false; + } + } + } /// /// The correlation matrix used for modeling dependency between the marginal distributions. @@ -339,6 +351,7 @@ public void SetParameters(UnivariateDistributionBase[] distributions) { if (distributions == null) throw new ArgumentNullException(nameof(Distributions)); _distributions = distributions; + _parametersValid = ValidateParameters(Array.Empty(), false) is null; _momentsComputed = false; _empiricalCDFCreated = false; _mvnCreated = false; @@ -356,15 +369,25 @@ public void SetParameters(IUnivariateDistribution[] distributions) { _distributions[i] = (UnivariateDistributionBase)distributions[i]; } + _parametersValid = ValidateParameters(Array.Empty(), false) is null; _momentsComputed = false; _empiricalCDFCreated = false; _mvnCreated = false; } /// + /// Thrown when the flattened parameter count does not match the component distributions. public override void SetParameters(IList parameters) { - if (Distributions == null || Distributions.Count == 0) return; + if (Distributions == null || Distributions.Count == 0) + { + _parametersValid = false; + return; + } + if (parameters.Count != NumberOfParameters) + { + throw new ArgumentException("The length of the parameter array is invalid.", nameof(parameters)); + } int t = 0; for (int i = 0; i < Distributions.Count; i++) @@ -378,6 +401,7 @@ public override void SetParameters(IList parameters) t += Distributions[i].NumberOfParameters; } + _parametersValid = ValidateParameters(parameters, false) is null; _momentsComputed = false; _empiricalCDFCreated = false; _mvnCreated = false; @@ -386,7 +410,12 @@ public override void SetParameters(IList parameters) /// public override ArgumentOutOfRangeException? ValidateParameters(IList parameters, bool throwException) { - if (Distributions.Count == 0) return new ArgumentOutOfRangeException(nameof(Distributions), "There must be at least 1 distribution"); + if (Distributions.Count == 0) + { + var exception = new ArgumentOutOfRangeException(nameof(Distributions), "There must be at least 1 distribution."); + if (throwException) throw exception; + return exception; + } for (int i = 0; i < Distributions.Count; i++) { if (Distributions[i].ParametersValid == false) @@ -1022,7 +1051,6 @@ private void CreateMultivariateNormal() var sigma = new double[D, D]; if (Dependency == Probability.DependencyType.PerfectlyNegative) { - CorrelationMatrix = new double[D, D]; double rho = -1d / (D - 1d) + Math.Sqrt(Tools.DoubleMachineEpsilon); for (int i = 0; i < D; i++) { diff --git a/Numerics/Distributions/Univariate/Deterministic.cs b/Numerics/Distributions/Univariate/Deterministic.cs index dcdbd9ef..75966e6c 100644 --- a/Numerics/Distributions/Univariate/Deterministic.cs +++ b/Numerics/Distributions/Univariate/Deterministic.cs @@ -190,6 +190,7 @@ public void Estimate(IList sample, ParameterEstimationMethod estimationM public void SetParameters(double value) { Value = value; + _parametersValid = ValidateParameters([value], false) is null; } /// diff --git a/Numerics/Distributions/Univariate/EmpiricalDistribution.cs b/Numerics/Distributions/Univariate/EmpiricalDistribution.cs index abfd33b9..4dc5f2bd 100644 --- a/Numerics/Distributions/Univariate/EmpiricalDistribution.cs +++ b/Numerics/Distributions/Univariate/EmpiricalDistribution.cs @@ -24,8 +24,8 @@ namespace Numerics.Distributions /// /// This distribution specifies a cumulative distribution with n points. The range of the distribution /// is set by the minimum and maximum arguments. Each point on the cumulative curve has an X value and - /// a probability. Points on the cumulative curve must be entered with increasing value and increasing - /// probability. Even though the (X,p) points define the distribution, and value between the minimum + /// a probability. Points on the cumulative curve must be entered with nondecreasing values and increasing + /// probabilities. Even though the (X,p) points define the distribution, any value between the minimum /// and maximum can be returned. /// /// @@ -102,7 +102,7 @@ public EmpiricalDistribution(OrderedPairedData orderedPairedData) } if (isAsc == true) { - opd = new OrderedPairedData(_xValues, _pValues, true, SortOrder.Ascending, true, SortOrder.Ascending); + opd = new OrderedPairedData(_xValues, _pValues, orderedPairedData.StrictX, SortOrder.Ascending, true, SortOrder.Ascending); } else { @@ -113,20 +113,22 @@ public EmpiricalDistribution(OrderedPairedData orderedPairedData) { _pValues = orderedPairedData.Select(v => v.Y).ToArray(); } + _parametersValid = ValidateData(opd, _pValues, false) is null; _momentsComputed = false; } /// /// Constructs a Univariate Empirical CDF from sample data. /// - /// + /// The sample values used as empirical ordinates; duplicate values are permitted. /// The plotting position formula type. Default = Weibull. public EmpiricalDistribution(IList sample, PlottingPositions.PlottingPostionType plottingPostionType = PlottingPositions.PlottingPostionType.Weibull) { _xValues = sample.ToArray(); Array.Sort(_xValues); _pValues = PlottingPositions.Function(_xValues.Count(), plottingPostionType)!; - opd = new OrderedPairedData(_xValues, _pValues, true, SortOrder.Ascending, true, SortOrder.Ascending); + opd = new OrderedPairedData(_xValues, _pValues, false, SortOrder.Ascending, true, SortOrder.Ascending); + _parametersValid = ValidateData(opd, _pValues, false) is null; _momentsComputed = false; } @@ -138,7 +140,7 @@ public EmpiricalDistribution(IList sample, PlottingPositions.PlottingPos /// /// Returns the array of X values. Points On the cumulative curve are specified - /// with increasing value and increasing probability. + /// with nondecreasing values and increasing probabilities. /// public ReadOnlyCollection XValues => new(_xValues); @@ -345,11 +347,17 @@ public IUnivariateDistribution Bootstrap(ParameterEstimationMethod estimationMet /// /// Array of X values. /// Array of probability values. Range 0 ≤ p ≤ 1. + /// The value and probability collections have different lengths. public void SetParameters(IList xValues, IList pValues) { + if (xValues.Count != pValues.Count) + { + throw new ArgumentException("The value and probability arrays must have the same length.", nameof(pValues)); + } _xValues = xValues.ToArray(); _pValues = pValues.ToArray(); - opd = new OrderedPairedData(xValues, pValues, true, SortOrder.Ascending, true, SortOrder.Ascending); + opd = new OrderedPairedData(xValues, pValues, false, SortOrder.Ascending, true, SortOrder.Ascending); + _parametersValid = ValidateData(opd, pValues, false) is null; _momentsComputed = false; } @@ -360,14 +368,54 @@ public void SetParameters(IList xValues, IList pValues) /// Array of probability values. Range 0 ≤ p ≤ 1. /// Sort order of X values. /// Sort order of probability values. + /// The value and probability collections have different lengths. public void SetParameters(IList xValues, IList pValues, SortOrder XOrder, SortOrder probabilityOrder) { + if (xValues.Count != pValues.Count) + { + throw new ArgumentException("The value and probability arrays must have the same length.", nameof(pValues)); + } _xValues = xValues.ToArray(); _pValues = pValues.ToArray(); - opd = new OrderedPairedData(xValues, pValues, true, XOrder, true, probabilityOrder); + opd = new OrderedPairedData(xValues, pValues, false, XOrder, true, probabilityOrder); + _parametersValid = ValidateData(opd, pValues, false) is null; _momentsComputed = false; } + /// + /// Validates empirical ordinates and their cumulative probabilities. + /// + /// The ordered paired data to validate. + /// The cumulative probabilities. + /// Whether to throw the validation exception. + /// when the data are valid; otherwise, the validation exception. + /// Thrown when is true and the data are invalid. + private static ArgumentOutOfRangeException? ValidateData(OrderedPairedData data, IList probabilities, bool throwException) + { + ArgumentOutOfRangeException? exception = null; + if (data.Count < 2) + { + exception = new ArgumentOutOfRangeException(nameof(data), "At least two empirical ordinates are required."); + } + else if (!data.IsValid) + { + exception = new ArgumentOutOfRangeException(nameof(data), "Empirical ordinates must be finite and follow the configured ordering."); + } + else + { + for (int i = 0; i < probabilities.Count; i++) + { + if (double.IsNaN(probabilities[i]) || double.IsInfinity(probabilities[i]) || probabilities[i] < 0d || probabilities[i] > 1d) + { + exception = new ArgumentOutOfRangeException(nameof(probabilities), "Empirical probabilities must be finite and between zero and one."); + break; + } + } + } + if (throwException && exception is not null) throw exception; + return exception; + } + /// public override void SetParameters(IList parameters) { @@ -427,6 +475,7 @@ public double PDF(double xl, double xu) /// public override double CDF(double X) { + if (_parametersValid == false) ValidateData(opd, _pValues, true); double p = 0; if (opd.OrderY == SortOrder.Ascending || opd.OrderY == SortOrder.None) { @@ -443,6 +492,7 @@ public override double CDF(double X) /// public override double InverseCDF(double probability) { + if (_parametersValid == false) ValidateData(opd, _pValues, true); // Validate probability if (probability < 0.0d || probability > 1.0d) throw new ArgumentOutOfRangeException(nameof(probability), "Probability must be between 0 and 1."); diff --git a/Numerics/Distributions/Univariate/GammaDistribution.cs b/Numerics/Distributions/Univariate/GammaDistribution.cs index 75d49291..196de9c7 100644 --- a/Numerics/Distributions/Univariate/GammaDistribution.cs +++ b/Numerics/Distributions/Univariate/GammaDistribution.cs @@ -689,12 +689,17 @@ public static double FrequencyFactorKp(double skewness, double probability) /// /// Coefficient of skewness. /// Probability between 0 and 1. + /// The partial derivative of the frequency factor with respect to skewness. public static double PartialKp(double skewness, double probability) { double C = skewness; double absC = Math.Abs(C); - // If skew is sufficiently close to zero, return standard Normal Z variate. - if (absC < 0.0001d) return Normal.StandardZ(probability); + // Use the Cornish-Fisher derivative limit at zero skew. + if (absC < 0.0001d) + { + double z = Normal.StandardZ(probability); + return (z * z - 1.0d) / 6.0d; + } // If abs(skew) is less than or equal to 2, use Cornish-Fisher transformation (Fisher and Cornish, 1960) if (absC <= 2d) diff --git a/Numerics/Distributions/Univariate/GeneralizedBeta.cs b/Numerics/Distributions/Univariate/GeneralizedBeta.cs index 33285244..81b236a4 100644 --- a/Numerics/Distributions/Univariate/GeneralizedBeta.cs +++ b/Numerics/Distributions/Univariate/GeneralizedBeta.cs @@ -241,7 +241,10 @@ public override double Mode { get { - if (Alpha <= 1.0d && Beta <= 1.0d) return (Min + Max) / 2.0d; + if (Alpha == 1.0d && Beta == 1.0d) return (Min + Max) / 2.0d; + if (Alpha < 1.0d && Beta < 1.0d) return (Min + Max) / 2.0d; + if (Alpha <= 1.0d && Beta >= 1.0d) return Min; + if (Alpha >= 1.0d && Beta <= 1.0d) return Max; double _mode = (Alpha - 1.0d) / (Alpha + Beta - 2.0d); return _mode * (Max - Min) + Min; } diff --git a/Numerics/Distributions/Univariate/GeneralizedLogistic.cs b/Numerics/Distributions/Univariate/GeneralizedLogistic.cs index 42eecc44..6c1087dc 100644 --- a/Numerics/Distributions/Univariate/GeneralizedLogistic.cs +++ b/Numerics/Distributions/Univariate/GeneralizedLogistic.cs @@ -465,8 +465,27 @@ public double[] ParametersFromLinearMoments(IList moments) double T3 = moments[2]; double T4 = moments[3]; double kappa = -T3; - double alpha = L2 * Math.Sin(kappa * Math.PI) / (kappa * Math.PI); - double xi = L1 - alpha * (1.0d / kappa - Math.PI / Math.Sin(kappa * Math.PI)); + double alpha; + double xi; + if (kappa == 0.0d) + { + alpha = L2; + xi = L1; + } + else if (Math.Abs(kappa) <= NearZero) + { + double kappa2 = kappa * kappa; + double pi2 = Math.PI * Math.PI; + double sinc = 1.0d - pi2 * kappa2 / 6.0d + pi2 * pi2 * kappa2 * kappa2 / 120.0d; + double reciprocalDifference = -pi2 * kappa / 6.0d - 7.0d * pi2 * pi2 * kappa * kappa2 / 360.0d; + alpha = L2 * sinc; + xi = L1 - alpha * reciprocalDifference; + } + else + { + alpha = L2 * Math.Sin(kappa * Math.PI) / (kappa * Math.PI); + xi = L1 - alpha * (1.0d / kappa - Math.PI / Math.Sin(kappa * Math.PI)); + } return [xi, alpha, kappa]; } @@ -478,8 +497,27 @@ public double[] LinearMomentsFromParameters(IList parameters) double kappa = parameters[2]; if (Math.Abs(kappa) >= 1.0d) throw new ArgumentOutOfRangeException(nameof(Kappa), "L-moments can only be defined for -1 < kappa < 1."); - double L1 = xi + alpha * (1.0d / kappa - Math.PI / Math.Sin(kappa * Math.PI)); - double L2 = alpha * kappa * Math.PI / Math.Sin(kappa * Math.PI); + double L1; + double L2; + if (kappa == 0.0d) + { + L1 = xi; + L2 = alpha; + } + else if (Math.Abs(kappa) <= NearZero) + { + double kappa2 = kappa * kappa; + double pi2 = Math.PI * Math.PI; + double reciprocalDifference = -pi2 * kappa / 6.0d - 7.0d * pi2 * pi2 * kappa * kappa2 / 360.0d; + double reciprocalSinc = 1.0d + pi2 * kappa2 / 6.0d + 7.0d * pi2 * pi2 * kappa2 * kappa2 / 360.0d; + L1 = xi + alpha * reciprocalDifference; + L2 = alpha * reciprocalSinc; + } + else + { + L1 = xi + alpha * (1.0d / kappa - Math.PI / Math.Sin(kappa * Math.PI)); + L2 = alpha * kappa * Math.PI / Math.Sin(kappa * Math.PI); + } double T3 = -kappa; double T4 = (1.0d + 5.0d * Math.Pow(kappa, 2.0d)) / 6.0d; return [L1, L2, T3, T4]; diff --git a/Numerics/Distributions/Univariate/Gumbel.cs b/Numerics/Distributions/Univariate/Gumbel.cs index 08656a31..938d53c2 100644 --- a/Numerics/Distributions/Univariate/Gumbel.cs +++ b/Numerics/Distributions/Univariate/Gumbel.cs @@ -227,11 +227,9 @@ public IUnivariateDistribution Bootstrap(ParameterEstimationMethod estimationMet /// The scale parameter α (alpha). public void SetParameters(double location, double scale) { - // Validate parameters - _parametersValid = ValidateParameters(location, scale, false) is null; - // Set parameters - Xi = location; + _xi = location; _alpha = scale; + _parametersValid = ValidateParameters(location, scale, false) is null; } /// diff --git a/Numerics/Distributions/Univariate/InverseChiSquared.cs b/Numerics/Distributions/Univariate/InverseChiSquared.cs index b6905ff8..d972c17f 100644 --- a/Numerics/Distributions/Univariate/InverseChiSquared.cs +++ b/Numerics/Distributions/Univariate/InverseChiSquared.cs @@ -220,14 +220,15 @@ public override double[] MaximumOfParameters /// public override void SetParameters(IList parameters) { - DegreesOfFreedom = (int)parameters[0]; - Sigma = parameters[1]; + _degreesOfFreedom = (int)parameters[0]; + _sigma = parameters[1]; + _parametersValid = ValidateParameters(parameters, false) is null; } /// public override ArgumentOutOfRangeException? ValidateParameters(IList parameters, bool throwException) { - if (parameters[0] < 1.0d) + if (double.IsNaN(parameters[0]) || double.IsInfinity(parameters[0]) || parameters[0] < 1.0d) { if (throwException) throw new ArgumentOutOfRangeException(nameof(DegreesOfFreedom), "The degrees of freedom ν (nu) must greater than or equal to one."); diff --git a/Numerics/Distributions/Univariate/KernelDensity.cs b/Numerics/Distributions/Univariate/KernelDensity.cs index 85e13bec..fe577c55 100644 --- a/Numerics/Distributions/Univariate/KernelDensity.cs +++ b/Numerics/Distributions/Univariate/KernelDensity.cs @@ -552,7 +552,7 @@ public override void SetParameters(IList parameters) /// private ArgumentOutOfRangeException? ValidateParameters(double value, bool throwException) { - if (value <= 0d) + if (double.IsNaN(value) || double.IsInfinity(value) || value <= 0d) { if (throwException) throw new ArgumentOutOfRangeException(nameof(Bandwidth), "The bandwidth must be a positive number!"); diff --git a/Numerics/Distributions/Univariate/LogPearsonTypeIII.cs b/Numerics/Distributions/Univariate/LogPearsonTypeIII.cs index f60e2f1e..4af02d19 100644 --- a/Numerics/Distributions/Univariate/LogPearsonTypeIII.cs +++ b/Numerics/Distributions/Univariate/LogPearsonTypeIII.cs @@ -579,6 +579,10 @@ public double[] ParametersFromLinearMoments(IList moments) double L1 = moments[0]; double L2 = moments[1]; double T3 = moments[2]; + if (T3 == 0.0d) + { + return [L1, L2 * Math.Sqrt(Math.PI), 0.0d]; + } var alpha = default(double); double z; // The following approximation has relative accuracy better than 5x10-5 for all values of alpha. @@ -595,14 +599,16 @@ public double[] ParametersFromLinearMoments(IList moments) double mu = L1; double gamma = 2.0d * Math.Pow(alpha, -0.5d) * Math.Sign(T3); - double sigma = L2 * Math.Pow(Math.PI, 0.5d) * Math.Pow(alpha, 0.5d) * Mathematics.SpecialFunctions.Gamma.Function(alpha) / Mathematics.SpecialFunctions.Gamma.Function(alpha + 0.5d); + double sigma; if (alpha < 100d) { - sigma = L2 * Math.Pow(Math.PI, 0.5d) * Math.Pow(alpha, 0.5d) * Mathematics.SpecialFunctions.Gamma.Function(alpha) / Mathematics.SpecialFunctions.Gamma.Function(alpha + 0.5d); + sigma = L2 * Math.Sqrt(Math.PI) * Math.Sqrt(alpha) * Mathematics.SpecialFunctions.Gamma.Function(alpha) / Mathematics.SpecialFunctions.Gamma.Function(alpha + 0.5d); } else { - sigma = Math.Sqrt(Math.PI) * L2 / (1d - 1d / (8.0d * alpha) + 1d / (128.0d * Math.Pow(alpha, 2d))); + double inverseAlpha = 1.0d / alpha; + double correction = 1.0d - inverseAlpha / 8.0d + inverseAlpha * inverseAlpha / 128.0d; + sigma = Math.Sqrt(Math.PI) * L2 / correction; } return [mu, sigma, gamma]; @@ -614,11 +620,24 @@ public double[] LinearMomentsFromParameters(IList parameters) double mu = parameters[0]; double sigma = parameters[1]; double gamma = parameters[2]; - double xi = mu - 2.0d * sigma / gamma; + if (gamma == 0.0d) + { + return [mu, sigma / Math.Sqrt(Math.PI), 0.0d, 0.12260172d]; + } double alpha = 4.0d / Math.Pow(gamma, 2d); double beta = 0.5d * sigma * gamma; - double L1 = xi + alpha * beta; - double L2 = Math.Abs(Math.Pow(Math.PI, -0.5d) * beta * Mathematics.SpecialFunctions.Gamma.Function(alpha + 0.5d) / Mathematics.SpecialFunctions.Gamma.Function(alpha)); + double L1 = mu; + double L2; + if (alpha < 100.0d) + { + L2 = Math.Abs(beta * Mathematics.SpecialFunctions.Gamma.Function(alpha + 0.5d) / (Math.Sqrt(Math.PI) * Mathematics.SpecialFunctions.Gamma.Function(alpha))); + } + else + { + double inverseAlpha = 1.0d / alpha; + double correction = 1.0d - inverseAlpha / 8.0d + inverseAlpha * inverseAlpha / 128.0d; + L2 = sigma / Math.Sqrt(Math.PI) * correction; + } // The following approximations are accurate to 10-6. double A0 = 0.32573501d; double A1 = 0.1686915d; @@ -656,6 +675,7 @@ public double[] LinearMomentsFromParameters(IList parameters) T3 = (1d + E1 * alpha + E2 * Math.Pow(alpha, 2d) + E3 * Math.Pow(alpha, 3d)) / (1d + F1 * alpha + F2 * Math.Pow(alpha, 2d) + F3 * Math.Pow(alpha, 3d)); T4 = (1d + G1 * alpha + G2 * Math.Pow(alpha, 2d) + G3 * Math.Pow(alpha, 3d)) / (1d + H1 * alpha + H2 * Math.Pow(alpha, 2d) + H3 * Math.Pow(alpha, 3d)); } + T3 *= Math.Sign(gamma); return [L1, L2, T3, T4]; } diff --git a/Numerics/Distributions/Univariate/Logistic.cs b/Numerics/Distributions/Univariate/Logistic.cs index 9a5f8099..ef9d9030 100644 --- a/Numerics/Distributions/Univariate/Logistic.cs +++ b/Numerics/Distributions/Univariate/Logistic.cs @@ -223,11 +223,9 @@ public IUnivariateDistribution Bootstrap(ParameterEstimationMethod estimationMet /// The scale parameter α (alpha). public void SetParameters(double location, double scale) { - // Validate parameters - _parametersValid = ValidateParameters(new[] { location, scale }, false) is null; - // Set parameters - Xi = location; + _xi = location; _alpha = scale; + _parametersValid = ValidateParameters(new[] { location, scale }, false) is null; } /// diff --git a/Numerics/Distributions/Univariate/Mixture.cs b/Numerics/Distributions/Univariate/Mixture.cs index 564eef67..763d022a 100644 --- a/Numerics/Distributions/Univariate/Mixture.cs +++ b/Numerics/Distributions/Univariate/Mixture.cs @@ -50,6 +50,8 @@ public Mixture(double[] weights, IUnivariateDistribution[] distributions) } private double[] _weights = null!; + private bool _isZeroInflated; + private double _zeroWeight; private UnivariateDistributionBase[] _distributions = null!; private EmpiricalDistribution _empiricalCDF = null!; private bool _momentsComputed = false; @@ -67,14 +69,47 @@ public Mixture(double[] weights, IUnivariateDistribution[] distributions) public UnivariateDistributionBase[] Distributions => _distributions; /// - /// Determines whether to use assign a weight to all data points less than or equal to zero. + /// Gets or sets whether a separate probability weight is assigned to values less than or equal to zero. /// - public bool IsZeroInflated { get; set; } = false; + public bool IsZeroInflated + { + get { return _isZeroInflated; } + set + { + _isZeroInflated = value; + RefreshConfigurationState(); + } + } /// - /// The zero-value weight used if the mixture is zero-inflated. + /// Gets or sets the zero-value probability weight used when the mixture is zero-inflated. /// - public double ZeroWeight { get; set; } = 0.0; + public double ZeroWeight + { + get { return _zeroWeight; } + set + { + _zeroWeight = value; + RefreshConfigurationState(); + } + } + + /// + /// Refreshes validity and cached results after zero-inflation configuration changes. + /// + private void RefreshConfigurationState() + { + if (_weights is null || _distributions is null) + { + _parametersValid = false; + } + else + { + _parametersValid = ValidateParameters(GetParameters, false) is null; + } + _momentsComputed = false; + _empiricalCDFCreated = false; + } /// /// Determines the interpolation transform for the X-values. @@ -377,6 +412,7 @@ public void SetParameters(double[] weights, UnivariateDistributionBase[] distrib _weights = weights; _distributions = distributions; + _parametersValid = ValidateParameters(GetParameters, false) is null; _momentsComputed = false; _empiricalCDFCreated = false; } @@ -399,6 +435,7 @@ public void SetParameters(double[] weights, IUnivariateDistribution[] distributi { _distributions[i] = (UnivariateDistributionBase)distributions[i]; } + _parametersValid = ValidateParameters(GetParameters, false) is null; _momentsComputed = false; _empiricalCDFCreated = false; } @@ -432,6 +469,7 @@ public void SetParameters(double[] weights, double[] parameters) Distributions[i].SetParameters(parms); t += Distributions[i].NumberOfParameters; } + _parametersValid = ValidateParameters(GetParameters, false) is null; _momentsComputed = false; _empiricalCDFCreated = false; } @@ -712,9 +750,27 @@ double[] MStep(double[] x) { wgt += likelihood[i, k]; } - } + } mleWeights[k] = wgt / N; } + + // Keep component weights on the configured simplex. For a zero-inflated + // mixture, finite-sample zero counts need not equal the fixed ZeroWeight. + double componentWeightSum = mleWeights.Sum(); + double componentWeightTarget = IsZeroInflated ? 1d - ZeroWeight : 1d; + if (componentWeightSum > 0d) + { + double scale = componentWeightTarget / componentWeightSum; + for (int k = 0; k < K; k++) + { + mleWeights[k] *= scale; + } + } + else + { + double weight = componentWeightTarget / K; + for (int k = 0; k < K; k++) mleWeights[k] = weight; + } // MLE var solver = new NelderMead(logLH, Np, x, Lowers, Uppers); solver.Maximize(); diff --git a/Numerics/Distributions/Univariate/NoncentralT.cs b/Numerics/Distributions/Univariate/NoncentralT.cs index 2b7c176e..0981a6e8 100644 --- a/Numerics/Distributions/Univariate/NoncentralT.cs +++ b/Numerics/Distributions/Univariate/NoncentralT.cs @@ -230,8 +230,9 @@ public override double[] MaximumOfParameters /// The noncentrality parameter μ (mu). public void SetParameters(double v, double mu) { - DegreesOfFreedom = v; - Noncentrality = mu; + _degreesOfFreedom = v; + _noncentrality = mu; + _parametersValid = ValidateParameters(v, mu, false) is null; } /// @@ -245,10 +246,12 @@ public override void SetParameters(IList parameters) /// /// The degrees of freedom ν (nu). Range: ν > 0. /// The noncentrality parameter μ (mu). - /// + /// Whether to throw the validation exception. + /// when the parameters are valid; otherwise, the validation exception. + /// Thrown when is true and either parameter is invalid. public ArgumentOutOfRangeException? ValidateParameters(double v, double mu, bool throwException) { - if (v < 1.0d) + if (double.IsNaN(v) || double.IsInfinity(v) || v < 1.0d) { if (throwException) throw new ArgumentOutOfRangeException(nameof(DegreesOfFreedom), "The degrees of freedom ν (nu) must greater than or equal to one."); @@ -279,8 +282,8 @@ public override double PDF(double x) double v = DegreesOfFreedom; if (x != 0) { - double A = NCT_CDF(x * Math.Sqrt(1 + 2 / v), v + 2, u); - double B = NCT_CDF(x, v, u); + double A = EvaluateCdfWithFallback(x * Math.Sqrt(1 + 2 / v), v + 2, u); + double B = EvaluateCdfWithFallback(x, v, u); double C = v / x; return C * (A - B); } @@ -299,7 +302,7 @@ public override double CDF(double x) // Validate parameters if (_parametersValid == false) ValidateParameters(_degreesOfFreedom, Noncentrality, true); - return NCTDist(x, DegreesOfFreedom, Noncentrality); + return EvaluateCdfSeries(x, DegreesOfFreedom, Noncentrality); } /// @@ -307,7 +310,7 @@ public override double InverseCDF(double probability) { // Validate probability if (probability < 0.0d || probability > 1.0d) - throw new ArgumentOutOfRangeException("probability", "Probability must be between 0 and 1."); + throw new ArgumentOutOfRangeException(nameof(probability), "Probability must be between 0 and 1."); if (probability == 0.0d) return Minimum; if (probability == 1.0d) @@ -316,7 +319,7 @@ public override double InverseCDF(double probability) if (_parametersValid == false) ValidateParameters(_degreesOfFreedom, Noncentrality, true); // - return NCT_INV(probability, DegreesOfFreedom, Noncentrality); + return InverseCdf(probability, DegreesOfFreedom, Noncentrality); } @@ -328,22 +331,23 @@ public override double InverseCDF(double probability) /// A single point in the distribution range. /// The degrees of freedom. /// The noncentrality parameter. - private double NCT_CDF(double t, double df, double delta) + /// The cumulative probability, using a normal approximation if the series does not converge. + private static double EvaluateCdfWithFallback(double t, double df, double delta) { - double ANS; + double result; try { - ANS = NCTDist(t, df, delta); + result = EvaluateCdfSeries(t, df, delta); } catch (ArgumentException) { // If the series fails to converge, use normal approximation - double Z = (t * (1.0d - 1.0d / (4.0d * df)) - delta) / Math.Sqrt(1.0d + t * t / (2.0d * df)); - ANS = Normal.StandardCDF(Z); - if (ANS < 0d) ANS = 0d; - if (ANS > 1d) ANS = 1d; + double z = (t * (1.0d - 1.0d / (4.0d * df)) - delta) / Math.Sqrt(1.0d + t * t / (2.0d * df)); + result = Normal.StandardCDF(z); + if (result < 0d) result = 0d; + if (result > 1d) result = 1d; } - return ANS; + return result; } /// @@ -354,132 +358,98 @@ private double NCT_CDF(double t, double df, double delta) /// A single point in the distribution range. /// The degrees of freedom. /// The noncentrality parameter. + /// The cumulative probability at . + /// Thrown when the AS 243 series exceeds its iteration limit. /// /// The function is based on ALGORITHM AS 243 APPL. STATIST. (1989), VOL.38, NO. 1. /// Original FORTRAN code can be found at: /// http://people.sc.fsu.edu/~jburkardt/f77_src/asa243/asa243.html /// - private double NCTDist(double t, double df, double delta) + private static double EvaluateCdfSeries(double t, double df, double delta) { - - // REAL FUNCTION TNC(T, DF, DELTA, IFAULT) - // - // ALGORITHM AS 243 APPL. STATIST. (1989), VOL.38, NO. 1 - // - // Cumulative probability at T of the non-central t-distribution - // with DF degrees of freedom (may be fractional) and non-centrality - // parameter DELTA. - // - // Note - requires the following auxiliary routines - // ALOGAM (X) - ACM 291 or AS 245 - // BETAIN (X, A, B, ALBETA, IFAULT) - AS 63 (updated in ASR 19) - // ALNORM (X, UPPER) - AS 66 - // - // Translated by William A. Huber. www.quantdec.com - // double a; - double ALBETA; + double logBeta; double b; - double DEL; - var N = default(int); - double ERRBD; - double GEVEN; - double GODD; - double LAMBDA; - double P; - double q; - double RXB; - double s; - double TT; + double adjustedDelta; + int iteration = 0; + double errorBound; + double gEven; + double gOdd; + double lambda; + double pWeight; + double qWeight; + double oneMinusXPowerB; + double remainingWeight; double x; - double XEVEN; - double XODD; - double TNC; - bool NEGDEL; - // - // Note - ITRMAX and ERRMAX may be changed to suit one's needs. - // - const int ITRMAX = 10000; - const double Errmax = 0.000000001d; - - // DATA ITRMAX/100.1/, ERRMAX/1.E-06/ - // - // Constants - R2PI = 1/ {GAMMA(1.5) * SQRT(2)} = SQRT(2 / PI) - // ALNRPI = Ln(SQRT(Pi)) - // + double xEven; + double xOdd; + double result; + bool reflected; + const int maxIterations = 10000; + const double maxError = 0.000000001d; const double zero = 0d; const double half = 0.5d; const double one = 1.0d; const double two = 2.0d; - const double r2pi = 0.797884560802865d; - const double alnrpi = 0.5723649429247d; - TNC = zero; - TT = t; - DEL = delta; - NEGDEL = false; + const double sqrtTwoOverPi = 0.797884560802865d; + const double logSqrtPi = 0.5723649429247d; + result = zero; + adjustedDelta = delta; + reflected = false; if (t < zero) { - NEGDEL = true; - TT = -TT; - DEL = -DEL; + reflected = true; + adjustedDelta = -adjustedDelta; } - // - // Initialize twin series (Guenther, J. Statist. Computn. Simuln. - // vol.6, 199, 1978). - // + + // Initialize the twin series of Guenther (1978). x = t * t / (t * t + df); - if (x <= zero) - goto Twenty; - LAMBDA = DEL * DEL; - P = half * Math.Exp(-half * LAMBDA); - q = r2pi * P * DEL; - s = half - P; - a = half; - b = half * df; - RXB = Math.Pow(one - x, b); - ALBETA = alnrpi + Gamma.LogGamma(b) - Gamma.LogGamma(a + b); - XODD = Beta.IncompleteRatio(x, a, b, ALBETA); - GODD = two * RXB * Math.Exp(a * Math.Log(x) - ALBETA); - XEVEN = one - RXB; - GEVEN = b * x * RXB; - TNC = P * XODD + q * XEVEN; - // - // Repeat until convergence - // - N = 1; - do + if (x > zero) { - a = a + one; - XODD = XODD - GODD; - XEVEN = XEVEN - GEVEN; - GODD = GODD * x * (a + b - one) / a; - GEVEN = GEVEN * x * (a + b - half) / (a + half); - P = P * LAMBDA / (two * N); - q = q * LAMBDA / (two * N + one); - s = s - P; - N = N + 1; - TNC = TNC + P * XODD + q * XEVEN; - ERRBD = two * s * (XODD - GODD); + lambda = adjustedDelta * adjustedDelta; + pWeight = half * Math.Exp(-half * lambda); + qWeight = sqrtTwoOverPi * pWeight * adjustedDelta; + remainingWeight = half - pWeight; + a = half; + b = half * df; + oneMinusXPowerB = Math.Pow(one - x, b); + logBeta = logSqrtPi + Gamma.LogGamma(b) - Gamma.LogGamma(a + b); + xOdd = Beta.IncompleteRatio(x, a, b, logBeta); + gOdd = two * oneMinusXPowerB * Math.Exp(a * Math.Log(x) - logBeta); + xEven = one - oneMinusXPowerB; + gEven = b * x * oneMinusXPowerB; + result = pWeight * xOdd + qWeight * xEven; + iteration = 1; + do + { + a = a + one; + xOdd = xOdd - gOdd; + xEven = xEven - gEven; + gOdd = gOdd * x * (a + b - one) / a; + gEven = gEven * x * (a + b - half) / (a + half); + pWeight = pWeight * lambda / (two * iteration); + qWeight = qWeight * lambda / (two * iteration + one); + remainingWeight = remainingWeight - pWeight; + iteration = iteration + 1; + result = result + pWeight * xOdd + qWeight * xEven; + errorBound = two * remainingWeight * (xOdd - gOdd); + } + while (errorBound > maxError && iteration <= maxIterations); } - while (ERRBD > Errmax && N <= ITRMAX); - // - Twenty: - ; - if (N > ITRMAX) + if (iteration > maxIterations) { throw new ArgumentException("Max number of iterations were exceeded."); } - if (NEGDEL) + if (reflected) { - // TNC = one - TNC - TNC = Normal.StandardCDF(DEL) - TNC; + result = Normal.StandardCDF(adjustedDelta) - result; } else { - TNC = TNC + (1.0d - Normal.StandardCDF(DEL)); - } // Upper tail area of N(0,1) - return TNC; + result = result + (1.0d - Normal.StandardCDF(adjustedDelta)); + } + return result; } /// @@ -488,84 +458,86 @@ private double NCTDist(double t, double df, double delta) /// Probability between 0 and 1. /// The degrees of freedom. /// The noncentrality parameter. - private double NCT_INV(double p, double df, double delta) - { - double t0, t1, t2, y0, y1, tInc, Slope; - int iter; - //double tEstimate; - const double ytol = 0.0000001d; // Y-tolerance - const double xtol = 0.0000001d; // X-tolerance - //const double w = 1.5d; // Over-shooting parameter (>= 1) - const int iterMax = 50; - t0 = NCTInv0(p, df, delta); - //tEstimate = t0; // Default if we run into trouble - y0 = NCTDist(t0, df, delta) - p; - if (y0 > ytol) + /// The quantile associated with . + private static double InverseCdf(double p, double df, double delta) + { + double lower; + double upper; + double candidate; + double lowerError; + double upperError; + double increment; + double slope; + int iteration; + const double probabilityTolerance = 0.0000001d; + const double quantileTolerance = 0.0000001d; + const int maxIterations = 50; + lower = InitialInverseCdfEstimate(p, df, delta); + lowerError = EvaluateCdfSeries(lower, df, delta) - p; + if (lowerError > probabilityTolerance) { - tInc = -1.0d; + increment = -1.0d; } - else if (y0 < -ytol) + else if (lowerError < -probabilityTolerance) { - tInc = 1.0d; + increment = 1.0d; } else { - return t0; + return lower; } - // - // Find a bracket through overshooting. - // - t1 = t0 + tInc; - y1 = NCTDist(t1, df, delta) - p; - iter = 0; - while ((y0 < 0d) != (y1 > 0d) && Math.Abs(t1 - t0) > xtol && iter < iterMax) + + // Find a bracket by secant extrapolation with a factor-of-two overshoot. + upper = lower + increment; + upperError = EvaluateCdfSeries(upper, df, delta) - p; + iteration = 0; + while ((lowerError < 0d) != (upperError > 0d) && Math.Abs(upper - lower) > quantileTolerance && iteration < maxIterations) { - // - // Use secant method to extrapolate a zero, but overshoot by w >= 1. - // - Slope = (y1 - y0) / (t1 - t0); - if (Slope == 0d) + slope = (upperError - lowerError) / (upper - lower); + if (slope == 0d) { - return (t1 + t0) / 2.0d; + return (upper + lower) / 2.0d; } - t2 = t0 - 2.0d * y0 / Slope; - t0 = t1; - t1 = t2; - y0 = y1; - y1 = NCTDist(t1, df, delta) - p; - iter = iter + 1; + candidate = lower - 2.0d * lowerError / slope; + lower = upper; + upper = candidate; + lowerError = upperError; + upperError = EvaluateCdfSeries(upper, df, delta) - p; + iteration = iteration + 1; } - // Solve for T using Brent - double ANS = Brent.Solve(x => NCTDist(x, df, delta) - p, Math.Min(t0, t1), Math.Max(t0, t1), xtol, reportFailure: false); - return ANS; + return Brent.Solve(x => EvaluateCdfSeries(x, df, delta) - p, Math.Min(lower, upper), Math.Max(lower, upper), quantileTolerance, reportFailure: false); } - private double NCTInv0(double P, double N, double D) + /// + /// Computes the Johnson and Kotz starting estimate for noncentral-t inversion. + /// + /// The cumulative probability. + /// The degrees of freedom. + /// The noncentrality parameter. + /// An initial quantile estimate for the root search. + /// Uses the Jennett and Welch approximation and falls back to an offset central-t quantile when necessary. + private static double InitialInverseCdfEstimate(double probability, double degreesOfFreedom, double noncentrality) { // Approximates percentage points of the non-central t distribution. - // P is the percentage, N is the degrees of freedom, D is the non-centrality parameter. // Source: Johnson & Kotz, Continuous Univariate Distributions, Volume 2. - double z = Normal.StandardZ(P); - // + double z = Normal.StandardZ(probability); + // Jennett & Welch approximation, formula (14.1). - // Intended for large values of D^2, such as are used for most tolerance interval calculations. - // - double b = Math.Exp(Gamma.LogGamma((N + 1d) / 2d) - Gamma.LogGamma(N / 2d)) * Math.Sqrt(2d / N); - double u2 = z * z; - double b2 = b * b; - double denom = b2 - u2 * (1d - b2); - double disc = b2 + (1d - b2) * (D * D - u2); - if (disc > 0d && Math.Abs(denom) > 1e-12) + // Intended for the large noncentralities common in tolerance-interval calculations. + double gammaRatio = Math.Exp(Gamma.LogGamma((degreesOfFreedom + 1d) / 2d) - Gamma.LogGamma(degreesOfFreedom / 2d)) * Math.Sqrt(2d / degreesOfFreedom); + double zSquared = z * z; + double gammaRatioSquared = gammaRatio * gammaRatio; + double denominator = gammaRatioSquared - zSquared * (1d - gammaRatioSquared); + double discriminant = gammaRatioSquared + (1d - gammaRatioSquared) * (noncentrality * noncentrality - zSquared); + if (discriminant > 0d && Math.Abs(denominator) > 1e-12) { - return (D * b + z * Math.Sqrt(disc)) / denom; - } - else - { - // Fallback: offset central t quantile by the noncentrality parameter - var st = new StudentT(N); - return st.InverseCDF(P) + D; + return (noncentrality * gammaRatio + z * Math.Sqrt(discriminant)) / denominator; } + + // Fall back to an offset central-t quantile. + var studentT = new StudentT(degreesOfFreedom); + return studentT.InverseCDF(probability) + noncentrality; } /// diff --git a/Numerics/Distributions/Univariate/PearsonTypeIII.cs b/Numerics/Distributions/Univariate/PearsonTypeIII.cs index 651c8dd7..6065f55b 100644 --- a/Numerics/Distributions/Univariate/PearsonTypeIII.cs +++ b/Numerics/Distributions/Univariate/PearsonTypeIII.cs @@ -409,6 +409,10 @@ public double[] ParametersFromLinearMoments(IList moments) double L1 = moments[0]; double L2 = moments[1]; double T3 = moments[2]; + if (T3 == 0.0d) + { + return [L1, L2 * Math.Sqrt(Math.PI), 0.0d]; + } double alpha = double.NaN; double z; // The following approximation has relative accuracy better than 5x10-5 for all values of alpha. @@ -424,7 +428,17 @@ public double[] ParametersFromLinearMoments(IList moments) } double gamma = 2.0d * Math.Pow(alpha, -0.5d) * Math.Sign(T3); - double sigma = L2 * Math.Pow(Math.PI, 0.5d) * Math.Pow(alpha, 0.5d) * Mathematics.SpecialFunctions.Gamma.Function(alpha) / Mathematics.SpecialFunctions.Gamma.Function(alpha + 0.5d); + double sigma; + if (alpha < 100.0d) + { + sigma = L2 * Math.Sqrt(Math.PI) * Math.Sqrt(alpha) * Mathematics.SpecialFunctions.Gamma.Function(alpha) / Mathematics.SpecialFunctions.Gamma.Function(alpha + 0.5d); + } + else + { + double inverseAlpha = 1.0d / alpha; + double correction = 1.0d - inverseAlpha / 8.0d + inverseAlpha * inverseAlpha / 128.0d; + sigma = Math.Sqrt(Math.PI) * L2 / correction; + } double mu = L1; return [mu, sigma, gamma]; } @@ -435,11 +449,24 @@ public double[] LinearMomentsFromParameters(IList parameters) double mu = parameters[0]; double sigma = parameters[1]; double gamma = parameters[2]; - double xi = mu - 2.0d * sigma / gamma; + if (gamma == 0.0d) + { + return [mu, sigma / Math.Sqrt(Math.PI), 0.0d, 0.12260172d]; + } double alpha = 4.0d / Math.Pow(gamma, 2d); double beta = 0.5d * sigma * gamma; - double L1 = xi + alpha * beta; - double L2 = Math.Abs(Math.Pow(Math.PI, -0.5d) * beta * Mathematics.SpecialFunctions.Gamma.Function(alpha + 0.5d) / Mathematics.SpecialFunctions.Gamma.Function(alpha)); + double L1 = mu; + double L2; + if (alpha < 100.0d) + { + L2 = Math.Abs(beta * Mathematics.SpecialFunctions.Gamma.Function(alpha + 0.5d) / (Math.Sqrt(Math.PI) * Mathematics.SpecialFunctions.Gamma.Function(alpha))); + } + else + { + double inverseAlpha = 1.0d / alpha; + double correction = 1.0d - inverseAlpha / 8.0d + inverseAlpha * inverseAlpha / 128.0d; + L2 = sigma / Math.Sqrt(Math.PI) * correction; + } // The following approximations are accurate to 10-6. double A0 = 0.32573501d; double A1 = 0.1686915d; @@ -477,6 +504,7 @@ public double[] LinearMomentsFromParameters(IList parameters) T3 = (1d + E1 * alpha + E2 * Math.Pow(alpha, 2d) + E3 * Math.Pow(alpha, 3d)) / (1d + F1 * alpha + F2 * Math.Pow(alpha, 2d) + F3 * Math.Pow(alpha, 3d)); T4 = (1d + G1 * alpha + G2 * Math.Pow(alpha, 2d) + G3 * Math.Pow(alpha, 3d)) / (1d + H1 * alpha + H2 * Math.Pow(alpha, 2d) + H3 * Math.Pow(alpha, 3d)); } + T3 *= Math.Sign(gamma); return [L1, L2, T3, T4]; } diff --git a/Numerics/Distributions/Univariate/StudentT.cs b/Numerics/Distributions/Univariate/StudentT.cs index 31b3cd13..4ca83e09 100644 --- a/Numerics/Distributions/Univariate/StudentT.cs +++ b/Numerics/Distributions/Univariate/StudentT.cs @@ -293,9 +293,10 @@ public override double[] MaximumOfParameters /// The degrees of freedom ν (nu). Range: ν > 0. public void SetParameters(double mu, double sigma, double v) { - Mu = mu; - Sigma = sigma; - DegreesOfFreedom = v; + _mu = mu; + _sigma = sigma; + _degreesOfFreedom = v; + _parametersValid = ValidateParameters(mu, sigma, v, false) is null; } /// @@ -325,7 +326,7 @@ public override void SetParameters(IList parameters) throw new ArgumentOutOfRangeException(nameof(Sigma), "Standard deviation must be positive."); return new ArgumentOutOfRangeException(nameof(Sigma), "Standard deviation must be positive."); } - if (v < 1.0d) + if (double.IsNaN(v) || double.IsInfinity(v) || v < 1.0d) { if (throwException) throw new ArgumentOutOfRangeException(nameof(DegreesOfFreedom), "The degrees of freedom ν (nu) must greater than or equal to one."); @@ -346,12 +347,13 @@ public override double PDF(double x) // Validate parameters if (_parametersValid == false) ValidateParameters(Mu, Sigma, DegreesOfFreedom, true); - double Z = (x - Mu) / Sigma; + double inverseSigma = 1.0d / Sigma; + double Z = (x - Mu) * inverseSigma; if (DegreesOfFreedom >= 100000000.0d) { - return Normal.StandardPDF(Z); + return Normal.StandardPDF(Z) * inverseSigma; } - return Math.Exp(Gamma.LogGamma((DegreesOfFreedom + 1.0d) / 2.0d) - Gamma.LogGamma(DegreesOfFreedom / 2.0d)) * Math.Pow(1.0d + Z * Z / DegreesOfFreedom, -0.5d * (DegreesOfFreedom + 1.0d)) / Math.Sqrt(DegreesOfFreedom * Math.PI); + return Math.Exp(Gamma.LogGamma((DegreesOfFreedom + 1.0d) / 2.0d) - Gamma.LogGamma(DegreesOfFreedom / 2.0d)) * Math.Pow(1.0d + Z * Z / DegreesOfFreedom, -0.5d * (DegreesOfFreedom + 1.0d)) * inverseSigma / Math.Sqrt(DegreesOfFreedom * Math.PI); } /// @@ -422,12 +424,7 @@ public override double InverseCDF(double probability) } double z = Beta.IncompleteInverse(0.5d * DegreesOfFreedom, 0.5d, 2.0d * probability); - double t = Math.Sqrt(DegreesOfFreedom / z - DegreesOfFreedom); - if (double.MaxValue * z < DegreesOfFreedom) - { - return rflg * double.MaxValue; - } - + double t = Math.Sqrt(DegreesOfFreedom) * Math.Sqrt(1.0d - z) / Math.Sqrt(z); return Mu + Sigma * (rflg * t); } } diff --git a/Numerics/Distributions/Univariate/TruncatedDistribution.cs b/Numerics/Distributions/Univariate/TruncatedDistribution.cs index f01b5e2f..8f9afaa3 100644 --- a/Numerics/Distributions/Univariate/TruncatedDistribution.cs +++ b/Numerics/Distributions/Univariate/TruncatedDistribution.cs @@ -27,14 +27,16 @@ public TruncatedDistribution(UnivariateDistributionBase basDistribution, double _baseDist = basDistribution; _min = min; _max = max; - _Fmin = _baseDist.CDF(_min); - _Fmax = _baseDist.CDF(_max); + var parameters = _baseDist.GetParameters.ToList(); + parameters.Add(min); + parameters.Add(max); + _parametersValid = ValidateParametersCore(parameters, false, out _Fmin, out _Fmax) is null; _momentsComputed = false; } private readonly UnivariateDistributionBase _baseDist; - private double _min = double.NegativeInfinity; - private double _max = double.PositiveInfinity; + private double _min; + private double _max; private double _Fmin, _Fmax; private bool _momentsComputed = false; private double[] u = [double.NaN, double.NaN, double.NaN, double.NaN]; @@ -230,32 +232,73 @@ public override double[] MaximumOfParameters } /// + /// Thrown when the flattened parameter count does not match the base distribution and two truncation bounds. public override void SetParameters(IList parameters) { - _baseDist.SetParameters(parameters.ToArray().Subset(parameters.Count - 2)); + if (parameters.Count != NumberOfParameters) + { + throw new ArgumentException("The length of the parameter array is invalid.", nameof(parameters)); + } + var validation = ValidateParametersCore(parameters, false, out double lowerProbability, out double upperProbability); + _baseDist.SetParameters(parameters.ToArray().Subset(0, parameters.Count - 3)); _min = parameters[parameters.Count - 2]; _max = parameters[parameters.Count - 1]; - _Fmin = _baseDist.CDF(_min); - _Fmax = _baseDist.CDF(_max); + _Fmin = lowerProbability; + _Fmax = upperProbability; + _parametersValid = validation is null; _momentsComputed = false; } /// public override ArgumentOutOfRangeException? ValidateParameters(IList parameters, bool throwException) { - if (_baseDist is not null) _baseDist.ValidateParameters(parameters.ToArray().Subset(0, parameters.Count - 2), throwException); - if (double.IsNaN(Min) || double.IsNaN(Max) || double.IsInfinity(Min) || double.IsInfinity(Max) || Min >= Max) + return ValidateParametersCore(parameters, throwException, out _, out _); + } + + /// + /// Validates a flattened base-distribution parameter vector and its truncation bounds. + /// + /// The base parameters followed by the lower and upper truncation bounds. + /// Whether to throw the first validation exception. + /// The base-distribution CDF at the lower bound when valid; otherwise . + /// The base-distribution CDF at the upper bound when valid; otherwise . + /// when all parameters are valid; otherwise, the validation exception. + /// Thrown when is true and validation fails. + private ArgumentOutOfRangeException? ValidateParametersCore(IList parameters, bool throwException, out double lowerProbability, out double upperProbability) + { + lowerProbability = double.NaN; + upperProbability = double.NaN; + if (parameters.Count != NumberOfParameters) + { + var exception = new ArgumentOutOfRangeException(nameof(parameters), "The flattened parameter count is invalid."); + if (throwException) throw exception; + return exception; + } + var baseParameters = parameters.ToArray().Subset(0, parameters.Count - 3); + var baseValidation = _baseDist.ValidateParameters(baseParameters, false); + if (baseValidation is not null) { - if (throwException) - throw new ArgumentOutOfRangeException(nameof(Min), "The min must be less than the max."); - return new ArgumentOutOfRangeException(nameof(Min), "The min must be less than the max."); + if (throwException) _baseDist.ValidateParameters(baseParameters, true); + return baseValidation; } - if (Math.Abs(_Fmin - _Fmax) < 1e-15) + double min = parameters[parameters.Count - 2]; + double max = parameters[parameters.Count - 1]; + if (double.IsNaN(min) || double.IsNaN(max) || double.IsInfinity(min) || double.IsInfinity(max) || min >= max) { - if (throwException) - throw new ArgumentOutOfRangeException(nameof(Min), "Truncation interval has zero probability mass."); - return new ArgumentOutOfRangeException(nameof(Min), "Truncation interval has zero probability mass."); - } + var exception = new ArgumentOutOfRangeException(nameof(Min), "The min must be finite and less than the max."); + if (throwException) throw exception; + return exception; + } + var candidate = _baseDist.Clone(); + candidate.SetParameters(baseParameters); + lowerProbability = candidate.CDF(min); + upperProbability = candidate.CDF(max); + if (Math.Abs(lowerProbability - upperProbability) < 1e-15) + { + var exception = new ArgumentOutOfRangeException(nameof(Min), "Truncation interval has zero probability mass."); + if (throwException) throw exception; + return exception; + } return null; } diff --git a/Numerics/Mathematics/Optimization/Local/BrentSearch.cs b/Numerics/Mathematics/Optimization/Local/BrentSearch.cs index 7e0f3079..f0226d86 100644 --- a/Numerics/Mathematics/Optimization/Local/BrentSearch.cs +++ b/Numerics/Mathematics/Optimization/Local/BrentSearch.cs @@ -157,45 +157,83 @@ protected override void Optimize() /// /// Bracket the objective function minimum. /// - /// Starting step size. Default = 1E-2. - /// Expansion factor. Default = 2. + /// The finite, nonzero starting step size. A negative value searches toward decreasing coordinates first. Default = 1E-2. + /// The finite geometric expansion factor, which must be greater than one. Default = 2. + /// is zero or non-finite, or is non-finite or not greater than one. + /// The bracket is not found within and is . + /// The search leaves the finite range of and is . + /// The objective function returns and is . + /// + /// The first three trial coordinates are equally spaced. After each unsuccessful downhill step, + /// the signed step is multiplied by until the middle value is no greater + /// than either endpoint. Failed searches leave the existing bounds unchanged. + /// public void Bracket(double s = 1E-2, double k = 2d) { + if (!Tools.IsFinite(s) || s == 0d) + throw new ArgumentOutOfRangeException(nameof(s), "The starting step size must be finite and nonzero."); + if (!Tools.IsFinite(k) || k <= 1d) + throw new ArgumentOutOfRangeException(nameof(k), "The expansion factor must be finite and greater than one."); + double a = LowerBound, b = a + s; + if (!Tools.IsFinite(a) || !Tools.IsFinite(b)) + { + UpdateStatus(OptimizationStatus.Failure, new ArithmeticException("The initial bracketing coordinates must be finite.")); + return; + } + double fa = ObjectiveFunction(new double[] { a }); double fb = ObjectiveFunction(new double[] { b }); - double c, fc, temp; + if (double.IsNaN(fa) || double.IsNaN(fb)) + { + UpdateStatus(OptimizationStatus.Failure, new InvalidOperationException("The objective function returned NaN while initializing the bracket.")); + return; + } + if (fb > fa) { - temp = a; + double temp = a; a = b; b = temp; - temp = fa; - fa = fb; - fb = temp; + fb = fa; s *= -1; } - while (true) + + for (int iteration = 0; iteration < MaxIterations; iteration++) { - c = b + s; - fc = ObjectiveFunction(new double[] { c }); - if (fc > fb) break; - a = b; + double c = b + s; + if (!Tools.IsFinite(c)) + { + UpdateStatus(OptimizationStatus.Failure, new ArithmeticException("The bracketing search exceeded the finite range of double-precision coordinates.")); + return; + } + + double fc = ObjectiveFunction(new double[] { c }); + if (double.IsNaN(fc)) + { + UpdateStatus(OptimizationStatus.Failure, new InvalidOperationException("The objective function returned NaN while expanding the bracket.")); + return; + } + + if (fc >= fb) + { + LowerBound = Math.Min(a, c); + UpperBound = Math.Max(a, c); + return; + } + + a = b; b = c; - fa = fb; fb = fc; - } - if (a < c) - { - LowerBound = a; - UpperBound = c; - } - else - { - LowerBound = c; - UpperBound = a; + s *= k; + if (!Tools.IsFinite(s)) + { + UpdateStatus(OptimizationStatus.Failure, new ArithmeticException("The geometric bracketing step exceeded the finite range of double precision.")); + return; + } } + UpdateStatus(OptimizationStatus.MaximumIterationsReached); } } diff --git a/Numerics/Numerics.csproj b/Numerics/Numerics.csproj index 23427c6c..32852b4b 100644 --- a/Numerics/Numerics.csproj +++ b/Numerics/Numerics.csproj @@ -6,7 +6,7 @@ True true embedded - + Numerics is a free and open-source .NET library developed by USACE-RMC, providing numerical methods, probability distributions, statistical analysis, optimization, machine learning, and Bayesian inference tools for quantitative risk assessment in water resources engineering. C. Haden Smith; Woodrow L. Fields; Julian Gonzalez; Sadie Niblett; Brennan Beam; Brian Skahill https://github.com/USACE-RMC/Numerics @@ -28,11 +28,9 @@ - - CS1587 - 2.1.2 - Version 2.1.2 improves time-series download reliability and statistics utilities. Time-series downloads now prefer IPv4 while preserving IPv6 fallback, avoid unrelated connectivity preflights on restricted networks, normalize duplicate USGS/GHCN timestamps deterministically, and fix Series.Clear()/RemoveAt behavior with large or duplicate-valued series. This release also adds RunningStatistics.Clone() and fixes Combine()/+ so combined statistics do not alias input instances when one operand is empty. - 2.1.2.0 + 2.1.3 + Version 2.1.3 hardens numerical edge cases and distribution validation. This release enforces XML documentation warnings, fixes audited interpolation/search, histogram, probability, MVNDST/COVSRT, Brent bracketing, and atomic double-addition behavior, improves copula clone and parameter-validity handling, and adds MultivariateNormal non-throwing covariance mutation plus marginal and conditional utilities with expanded regression coverage. + 2.1.3.0 @@ -72,7 +70,7 @@ - + @@ -107,5 +105,5 @@ Resources.Designer.cs - + diff --git a/Numerics/Sampling/MCMC/Base/MCMCSampler.cs b/Numerics/Sampling/MCMC/Base/MCMCSampler.cs index bbbb8465..52d1f370 100644 --- a/Numerics/Sampling/MCMC/Base/MCMCSampler.cs +++ b/Numerics/Sampling/MCMC/Base/MCMCSampler.cs @@ -359,6 +359,7 @@ protected virtual void InitializeCustomSettings() { } /// /// Initialize the Markov Chains. /// + /// The initialized parameter state for each Markov chain. protected virtual ParameterSet[] InitializeChains() { if (Initialize == InitializationType.UserDefined) @@ -384,7 +385,6 @@ protected virtual ParameterSet[] InitializeChains() // Use differential evolution to find a global optimum var lowerBounds = PriorDistributions.Select(x => x.Minimum).ToArray(); var upperBounds = PriorDistributions.Select(x => x.Maximum).ToArray(); - var inititals = lowerBounds.Add(upperBounds).Divide(2d); var DE = new DifferentialEvolution((x) => { return LogLikelihoodFunction(x); }, NumberOfParameters, lowerBounds, upperBounds); DE.ReportFailure = false; DE.Maximize(); @@ -394,7 +394,7 @@ protected virtual ParameterSet[] InitializeChains() { _mapSuccessful = true; // Get MAP - MAP = DE.BestParameterSet.Clone(); + MAP = new ParameterSet((double[])DE.BestParameterSet.Values.Clone(), -DE.BestParameterSet.Fitness); // Get Fisher Information Matrix if (DE.Hessian == null) throw new InvalidOperationException("Hessian matrix is not available."); var fisher = DE.Hessian * -1d; diff --git a/Numerics/Sampling/MCMC/NUTS.cs b/Numerics/Sampling/MCMC/NUTS.cs index 76be0966..263cdc32 100644 --- a/Numerics/Sampling/MCMC/NUTS.cs +++ b/Numerics/Sampling/MCMC/NUTS.cs @@ -396,7 +396,7 @@ private void LeapfrogInPlace(double[] theta, double[] momentum, double epsilon, double[] invMass = _inverseMassMatrix[chainIndex]; // Half-step momentum update - double[] grad = NumericalDerivative.Gradient((y) => SafeLogLikelihood(y), theta, _lowerBounds, _upperBounds); + double[] grad = GradientFunction(theta).Array; for (int j = 0; j < D; j++) momentum[j] += grad[j] * halfEps; @@ -411,7 +411,7 @@ private void LeapfrogInPlace(double[] theta, double[] momentum, double epsilon, } // Half-step momentum update - grad = NumericalDerivative.Gradient((y) => SafeLogLikelihood(y), theta, _lowerBounds, _upperBounds); + grad = GradientFunction(theta).Array; for (int j = 0; j < D; j++) momentum[j] += grad[j] * halfEps; } diff --git a/Numerics/Utilities/Tools.cs b/Numerics/Utilities/Tools.cs index 3455d4f0..0b6215d5 100644 --- a/Numerics/Utilities/Tools.cs +++ b/Numerics/Utilities/Tools.cs @@ -3,6 +3,7 @@ using System.Collections.Generic; using System.IO; using System.IO.Compression; +using System.Runtime.CompilerServices; using System.Threading; namespace Numerics @@ -182,9 +183,11 @@ public static double Pow(double a, int b) /// Returns the base 10 logarithm of a specified number. /// /// The number whose logarithm is to be found. + /// The guarded base-10 logarithm of . + [System.Runtime.CompilerServices.MethodImpl(System.Runtime.CompilerServices.MethodImplOptions.AggressiveInlining)] public static double Log10(double x) { - if (x < 1E-16 && Math.Sign(x) != -1) x = 1E-16; + if (x < 1E-16 && x >= 0d) x = 1E-16; return Math.Log10(x); } @@ -662,17 +665,57 @@ public static double Max(IList values, IList indicators, bool useCo /// /// Supporting function used to add doubles from an interlocked parallel loop. /// - /// The value to add to. + /// The accumulator to update atomically. /// The value to add. + /// The accumulator value produced by the atomic addition. + /// + /// All concurrent writes to must use this method or another + /// interlocked operation. The compare-and-swap operation can retry under contention. Concurrent + /// updates are not lost, but their order depends on thread scheduling, so floating-point reductions + /// using this method are not guaranteed to be bit-reproducible. + /// + [MethodImpl(MethodImplOptions.AggressiveInlining)] public static double ParallelAdd(ref double valueToAddTo, double valueToAdd) { - double newCurrentValue = valueToAddTo; + double currentValue = valueToAddTo; + double newValue = currentValue + valueToAdd; + double observedValue = Interlocked.CompareExchange(ref valueToAddTo, newValue, currentValue); + if (observedValue == currentValue && + (currentValue != 0d || + BitConverter.DoubleToInt64Bits(observedValue) == BitConverter.DoubleToInt64Bits(currentValue))) + { + return newValue; + } + if (double.IsNaN(observedValue) && double.IsNaN(currentValue)) return newValue; + return ParallelAddRetry(ref valueToAddTo, valueToAdd, observedValue); + } + + /// + /// Retries an atomic double addition after a failed compare-and-swap operation. + /// + /// The accumulator to update atomically. + /// The value to add. + /// The accumulator value observed by the failed operation. + /// The accumulator value produced by the atomic addition. + /// + /// The observed value seeds the lock-free retry loop so that another non-volatile read is unnecessary. + /// Signed zeros require a representation comparison because numerical equality does not distinguish + /// positive zero from negative zero, while compare-and-swap does. + /// + private static double ParallelAddRetry(ref double valueToAddTo, double valueToAdd, double observedValue) + { + double newCurrentValue = observedValue; while (true) { double currentValue = newCurrentValue; double newValue = currentValue + valueToAdd; newCurrentValue = Interlocked.CompareExchange(ref valueToAddTo, newValue, currentValue); - if (newCurrentValue == currentValue) return newValue; + if (newCurrentValue == currentValue && + (currentValue != 0d || + BitConverter.DoubleToInt64Bits(newCurrentValue) == BitConverter.DoubleToInt64Bits(currentValue))) + { + return newValue; + } if (double.IsNaN(newCurrentValue) && double.IsNaN(currentValue)) return newValue; } } diff --git a/Test_Numerics/Data/Interpolation/Test_Bilinear.cs b/Test_Numerics/Data/Interpolation/Test_Bilinear.cs index a605e6b4..9619cb43 100644 --- a/Test_Numerics/Data/Interpolation/Test_Bilinear.cs +++ b/Test_Numerics/Data/Interpolation/Test_Bilinear.cs @@ -304,6 +304,7 @@ public void Test_BilinearEdgeCases() x2 = 75; y = bilinear.Interpolate(x1, x2); Assert.AreEqual(859.405, y, 1E-6); + // Bottom x1 = 600; x2 = 225; @@ -323,5 +324,38 @@ public void Test_BilinearEdgeCases() Assert.AreEqual(929.65000, y, 1E-6); } - } + /// + /// Verifies guarded logarithmic transforms for zero and sub-floor coordinates and ordinates. + /// + [TestMethod] + public void Test_LogarithmicFloorMatchesLinearInterpolation() + { + var coordinates = new[] { 0d, 1E-15, 1d }; + var secondCoordinates = new[] { 0d, 1E-15, 1d }; + var values = new[,] + { + { 0d, 0d, 0d }, + { 1E-15, 1E-15, 1E-15 }, + { 1d, 1d, 1d } + }; + var bilinear = new Bilinear(coordinates, secondCoordinates, values) + { + X1Transform = Transform.Logarithmic, + X2Transform = Transform.Logarithmic, + YTransform = Transform.Logarithmic + }; + var linear = new Linear(coordinates, new[] { 0d, 1E-15, 1d }) + { + XTransform = Transform.Logarithmic, + YTransform = Transform.Logarithmic + }; + + double atZero = bilinear.Interpolate(0d, 0d); + double belowFloor = bilinear.Interpolate(5E-17, 5E-17); + Assert.IsFalse(double.IsNaN(atZero) || double.IsInfinity(atZero)); + Assert.IsFalse(double.IsNaN(belowFloor) || double.IsInfinity(belowFloor)); + Assert.AreEqual(linear.Interpolate(5E-17), belowFloor, 1E-28); + } + } + } diff --git a/Test_Numerics/Data/Interpolation/Test_Linear.cs b/Test_Numerics/Data/Interpolation/Test_Linear.cs index 7accecb5..5aa923fc 100644 --- a/Test_Numerics/Data/Interpolation/Test_Linear.cs +++ b/Test_Numerics/Data/Interpolation/Test_Linear.cs @@ -1,4 +1,4 @@ -using System; +using System; using Microsoft.VisualStudio.TestTools.UnitTesting; using Numerics.Data; @@ -238,7 +238,7 @@ public void Test_Z_R() public void Test_RevLin() { var XArray = new double[] { 50d, 100d, 150d, 200d, 250d }; - var YArray = new double[] { 100d, 200d, 300d, 400d, 500d }; + var YArray = new double[] { 100d, 200d, 300d, 400d, 500d }; Array.Reverse(XArray); Array.Reverse(YArray); var LI = new Linear(XArray, YArray, SortOrder.Descending); @@ -300,6 +300,9 @@ public void Test_Rev_Z() } // ??? + /// + /// Tests linear interpolation from list inputs. + /// [TestMethod] public void Test_Lin_List() { diff --git a/Test_Numerics/Data/Interpolation/Test_Search.cs b/Test_Numerics/Data/Interpolation/Test_Search.cs new file mode 100644 index 00000000..bc8b7ff9 --- /dev/null +++ b/Test_Numerics/Data/Interpolation/Test_Search.cs @@ -0,0 +1,100 @@ +using System.Collections.Generic; +using System.Linq; +using Microsoft.VisualStudio.TestTools.UnitTesting; +using Numerics.Data; + +namespace Data.Interpolation +{ + /// + /// Regression tests for interpolation search helpers. + /// + [TestClass] + public class Test_Search + { + /// + /// Gets the interior, endpoint, and out-of-range values exercised by every search overload. + /// + private static readonly double[] Queries = { -5d, 0d, 5d, 25d, 35d, 40d, 45d }; + + /// + /// Verifies that ascending searches retain their existing behavior for every overload. + /// + [TestMethod] + public void AscendingSearchesAgreeWithSequential() + { + VerifyAllOverloads(new[] { 0d, 10d, 20d, 30d, 40d }, SortOrder.Ascending); + } + + /// + /// Verifies descending interiors, endpoints, and out-of-range values for every overload. + /// + [TestMethod] + public void DescendingSearchesAgreeWithSequential() + { + VerifyAllOverloads(new[] { 40d, 30d, 20d, 10d, 0d }, SortOrder.Descending); + } + + /// + /// Verifies nonzero bisection starts and hunt guesses in both directions. + /// + [TestMethod] + public void NonzeroStartsReturnTheSequentialBracket() + { + VerifyStarts(new[] { 0d, 10d, 20d, 30d, 40d }, SortOrder.Ascending, 35d, 15d); + VerifyStarts(new[] { 40d, 30d, 20d, 10d, 0d }, SortOrder.Descending, 15d, 25d); + } + + /// + /// Verifies that bisection and hunt agree with sequential search for every supported data container. + /// + /// The ordered values to search. + /// The sort direction of . + private static void VerifyAllOverloads(double[] values, SortOrder order) + { + var paired = new OrderedPairedData(values, values, true, order, true, order); + IList ordinates = values.Select(value => new Ordinate(value, value)).ToList(); + + foreach (double query in Queries) + { + int expected = Search.Sequential(query, values, 0, order); + Assert.AreEqual(expected, Search.Bisection(query, values, 0, order)); + Assert.AreEqual(expected, Search.Hunt(query, values, 0, order)); + + expected = Search.Sequential(query, paired); + Assert.AreEqual(expected, Search.Bisection(query, paired)); + Assert.AreEqual(expected, Search.Hunt(query, paired)); + + expected = Search.Sequential(query, ordinates, 0, order); + Assert.AreEqual(expected, Search.Bisection(query, ordinates, 0, order)); + Assert.AreEqual(expected, Search.Hunt(query, ordinates, 0, order)); + } + } + + /// + /// Verifies bisection and hunt from nonzero starting indices for every supported data container. + /// + /// The ordered values to search. + /// The sort direction of . + /// The value used to verify bisection. + /// The value used to verify hunt. + private static void VerifyStarts( + double[] values, + SortOrder order, + double bisectionQuery, + double huntQuery) + { + var paired = new OrderedPairedData(values, values, true, order, true, order); + IList ordinates = values.Select(value => new Ordinate(value, value)).ToList(); + + int expected = Search.Sequential(bisectionQuery, values, 0, order); + Assert.AreEqual(expected, Search.Bisection(bisectionQuery, values, 2, order)); + Assert.AreEqual(expected, Search.Bisection(bisectionQuery, paired, 2)); + Assert.AreEqual(expected, Search.Bisection(bisectionQuery, ordinates, 2, order)); + + expected = Search.Sequential(huntQuery, values, 0, order); + Assert.AreEqual(expected, Search.Hunt(huntQuery, values, 3, order)); + Assert.AreEqual(expected, Search.Hunt(huntQuery, paired, 3)); + Assert.AreEqual(expected, Search.Hunt(huntQuery, ordinates, 3, order)); + } + } +} diff --git a/Test_Numerics/Data/Regression/Test_LinearRegression.cs b/Test_Numerics/Data/Regression/Test_LinearRegression.cs index 26330223..cc71bc75 100644 --- a/Test_Numerics/Data/Regression/Test_LinearRegression.cs +++ b/Test_Numerics/Data/Regression/Test_LinearRegression.cs @@ -1,4 +1,4 @@ -using System.Collections.Generic; +using System.Collections.Generic; using System.Diagnostics; // keep to be able to uncomment and run Debug.WriteLine() using Microsoft.VisualStudio.TestTools.UnitTesting; using Numerics.Data; @@ -7,7 +7,7 @@ namespace Data.Regression { /// - /// Unit tests for the Linear Regression class. The linear regression models are validated against R's "lm( )" method, while + /// Unit tests for the Linear Regression class. The linear regression models are validated against R's "lm( )" method, while /// the prediction methods are validated with the "predict.lm( )" function. /// /// @@ -25,7 +25,7 @@ namespace Data.Regression public class Test_LinearRegression { /// - /// Test simple linear regression. + /// Test simple linear regression. /// [TestMethod] public void Test_SimpleLinearRegression() @@ -62,7 +62,7 @@ public void Test_SimpleLinearRegression() Assert.AreEqual(ar2, true_ar2, 1E-4); Assert.AreEqual(df, true_df); - // Test summary output table + // Test summary output table var summary = LM.Summary(); for (int i = 0; i < summary.Count; i++) { @@ -87,7 +87,7 @@ public void Test_SimpleLinearRegression() } /// - /// Test multiple linear regression. + /// Test multiple linear regression. /// [TestMethod] public void Test_MultipleLinearRegression() @@ -121,14 +121,14 @@ public void Test_MultipleLinearRegression() Assert.AreEqual(par[i], true_par[i], 1E-4); Assert.AreEqual(sig[i], true_sig[i], 1E-4); } - + Assert.AreEqual(se, true_se, 1E-4); Assert.AreEqual(r2, true_r2, 1E-3); Assert.AreEqual(ar2, true_ar2, 1E-3); Assert.AreEqual(df, true_df); - // Test summary output table + // Test summary output table var summary = LM.Summary(); for (int i = 0; i < summary.Count; i++) { @@ -136,23 +136,23 @@ public void Test_MultipleLinearRegression() // Debug.WriteLine(summary[i]); } - /** The below summary table is shown for demonstration. + /* The below summary table is shown for demonstration. * This is the same summary that would be outputted from the R lm() method. - * + * * Model for predicting Consumption: * Parameters: * Estimate Std. Error t value Pr(>|t|) * Intercept 0.26729 0.03721 7.184 1.68E-011 *** - * Income 0.71448 0.04219 16.934 <1E-15 *** + * Income 0.71448 0.04219 16.934 <1E-15 *** * Production 0.04589 0.02588 1.773 0.0778 . - * Savings -0.04527 0.00278 -16.287 <1E-15 *** + * Savings -0.04527 0.00278 -16.287 <1E-15 *** * Unemployment -0.20477 0.10550 -1.941 0.0538 . * --- * Signif.codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 * * Residual standard error: 0.3286 on 182 degrees of freedom * Multiple R-squared: 0.7540, Adjusted R-squared: 0.7486 - * F-statistic: 139.5, on 4 and 182 DF, p-value: < 1E-15 + * F-statistic: 139.5, on 4 and 182 DF, p-value: < 1E-15 * * Residuals: * Min 1Q Median 3Q Max @@ -279,7 +279,7 @@ public void Test_PredictionIntervals() { 0.1024060, 2.110207, 1.10630631 }, }; var testY = LM.PredictionIntervals(testX); - + for(int i = 0; i < 4; ++i) { for (int j = 0; j < 3; j++) @@ -306,7 +306,7 @@ public void Test_PredictionIntervals_Multi() var LM = new LinearRegression(x, y); var tx = new double[,] - { + { { -1, -1, 1, -0.25 }, { -2, -2, 2, -0.15 }, { 1, 1, 3, 0.15 }, diff --git a/Test_Numerics/Data/Statistics/Test_Histogram.cs b/Test_Numerics/Data/Statistics/Test_Histogram.cs index 5da9f3e6..41890b93 100644 --- a/Test_Numerics/Data/Statistics/Test_Histogram.cs +++ b/Test_Numerics/Data/Statistics/Test_Histogram.cs @@ -192,5 +192,44 @@ public void Test_Indexing() Assert.AreEqual(i, index); } } + /// + /// Verifies adaptive range extension, exact boundaries, batches, and unchanged in-range behavior. + /// + [TestMethod] + public void Test_AddData_AdaptsEndpointBins() + { + var histogram = new Histogram(new[] { 0d, 10d }, 2); + int firstFrequency = histogram[0].Frequency; + int lastFrequency = histogram[histogram.NumberOfBins - 1].Frequency; + + histogram.AddData(2d); + Assert.AreEqual(0d, histogram.LowerBound); + Assert.AreEqual(10d, histogram.UpperBound); + Assert.AreEqual(firstFrequency + 1, histogram[0].Frequency); + + histogram.AddData(0d); + histogram.AddData(10d); + Assert.AreEqual(firstFrequency + 2, histogram[0].Frequency); + Assert.AreEqual(lastFrequency + 1, histogram[histogram.NumberOfBins - 1].Frequency); + + histogram.AddData(-5d); + histogram.AddData(-10d); + histogram.AddData(15d); + histogram.AddData(20d); + Assert.AreEqual(-10d, histogram.LowerBound); + Assert.AreEqual(20d, histogram.UpperBound); + Assert.AreEqual(-10d, histogram[0].LowerBound); + Assert.AreEqual(20d, histogram[histogram.NumberOfBins - 1].UpperBound); + + histogram.AddData(new[] { -12d, 22d, 5d }); + Assert.AreEqual(-12d, histogram.LowerBound); + Assert.AreEqual(22d, histogram.UpperBound); + Assert.AreEqual(-12d, histogram[0].LowerBound); + Assert.AreEqual(22d, histogram[histogram.NumberOfBins - 1].UpperBound); + Assert.AreEqual(12, histogram.DataCount); + Assert.AreEqual( + histogram.DataCount, + Enumerable.Range(0, histogram.NumberOfBins).Sum(index => histogram[index].Frequency)); + } } } diff --git a/Test_Numerics/Data/Statistics/Test_HypothesisTests.cs b/Test_Numerics/Data/Statistics/Test_HypothesisTests.cs index 6e3ed6dc..e1624d47 100644 --- a/Test_Numerics/Data/Statistics/Test_HypothesisTests.cs +++ b/Test_Numerics/Data/Statistics/Test_HypothesisTests.cs @@ -1,4 +1,4 @@ -using System; +using System; using Microsoft.VisualStudio.TestTools.UnitTesting; using Numerics.Data.Statistics; using Numerics.Distributions; @@ -19,7 +19,7 @@ namespace Data.Statistics /// References: /// /// - /// R Core Team (2013). R: A language and environment for statistical computing. R Foundation for Statistical Computing, + /// R Core Team (2013). R: A language and environment for statistical computing. R Foundation for Statistical Computing, /// Vienna, Austria. ISBN 3-900051-07-0, URL http://www.R-project.org/. /// /// @@ -137,7 +137,7 @@ public void Test_FtestModels() /// /// References: /// - /// + /// /// Lukasz Komsta (2005). moments: Moments, Cumulants, Skewness, Kurtosis and Related Tests. R package version 0.14.1, https://cran.r-project.org/web/packages/moments /// /// @@ -313,7 +313,7 @@ public void Test_MultipleGrubbsBeck() /// /// /// References: - /// Scrucca, L., Fop, M., Murphy, T. B., & Raftery, A. E. (2016). mclust 5: Clustering, Classification and Density Estimation Using Gaussian Finite Mixture Models. The R Journal, 8(1), 289–317. https://doi.org/10.32614/RJ-2016-021 + /// Scrucca, L., Fop, M., Murphy, T. B., & Raftery, A. E. (2016). mclust 5: Clustering, Classification and Density Estimation Using Gaussian Finite Mixture Models. The R Journal, 8(1), 289–317. https://doi.org/10.32614/RJ-2016-021 /// [TestMethod] public void Test_UnimodalityTest() @@ -327,8 +327,8 @@ public void Test_UnimodalityTest() double true_pval = 0.4142441; Assert.AreEqual(true_pval, pval, 1E-4); - var bimodalData = new double[] { 3.8, 3.9, 4.0, 3.9, 4.0, 4.1, 4.0, 3.8, 3.9, 4.0, - 4.1, 4.0, 3.9, 4.0, 4.0, 4.3, 4.4, 4.4, 4.5, 4.4, 4.3, + var bimodalData = new double[] { 3.8, 3.9, 4.0, 3.9, 4.0, 4.1, 4.0, 3.8, 3.9, 4.0, + 4.1, 4.0, 3.9, 4.0, 4.0, 4.3, 4.4, 4.4, 4.5, 4.4, 4.3, 4.4, 4.5, 4.4, 4.3, 4.4, 4.4, 4.5, 4.4, 4.3 }; pval = HypothesisTests.UnimodalityTest(bimodalData); diff --git a/Test_Numerics/Data/Statistics/Test_Probability.cs b/Test_Numerics/Data/Statistics/Test_Probability.cs index 57b1cb1e..8b95c3c0 100644 --- a/Test_Numerics/Data/Statistics/Test_Probability.cs +++ b/Test_Numerics/Data/Statistics/Test_Probability.cs @@ -108,6 +108,39 @@ public void Test_JointABCD_Independent_PCM() } + /// + /// Verifies HPCM remains finite for endpoint and subnormal marginal probabilities. + /// + [TestMethod] + public void Test_JointProbabilityHPCM_ExtremeProbabilitiesRemainFinite() + { + var probabilities = new[] { 0d, 1E-320, 0.5d }; + var indicators = new[] { 1, 1, 1 }; + var correlation = new double[,] + { + { 1d, 0.5d, 0.25d }, + { 0.5d, 1d, 0.25d }, + { 0.25d, 0.25d, 1d } + }; + var conditionalProbabilities = new double[probabilities.Length]; + + double joint = Probability.JointProbabilityHPCM( + probabilities, + indicators, + correlation, + conditionalProbabilities); + + Assert.IsFalse(double.IsNaN(joint)); + Assert.IsFalse(double.IsInfinity(joint)); + Assert.IsTrue(joint >= 0d && joint <= 1d); + foreach (double conditionalProbability in conditionalProbabilities) + { + Assert.IsFalse(double.IsNaN(conditionalProbability)); + Assert.IsFalse(double.IsInfinity(conditionalProbability)); + Assert.IsTrue(conditionalProbability >= 0d && conditionalProbability <= 1d); + } + } + /// /// Test joint probability of ABCD assuming perfect positive dependence. /// diff --git a/Test_Numerics/Distributions/Bivariate Copulas/Test_AMHCopula.cs b/Test_Numerics/Distributions/Bivariate Copulas/Test_AMHCopula.cs index ada6f7f4..9f06993f 100644 --- a/Test_Numerics/Distributions/Bivariate Copulas/Test_AMHCopula.cs +++ b/Test_Numerics/Distributions/Bivariate Copulas/Test_AMHCopula.cs @@ -176,6 +176,45 @@ public void Test_TailDependence() Assert.AreEqual(0.0, copulaNeg.UpperTailDependence, 1E-10, "AMH copula should have no upper tail dependence for negative θ."); Assert.AreEqual(0.0, copulaNeg.LowerTailDependence, 1E-10, "AMH copula should have no lower tail dependence for negative θ."); } + + /// + /// Test that ParametersValid tracks the dependency parameter's valid range. + /// + [TestMethod] + public void Test_ParametersValid() + { + var copula = new AMHCopula(0.4); + Assert.IsTrue(copula.ParametersValid); + + copula.Theta = 1.5; // above the maximum of 1 + Assert.IsFalse(copula.ParametersValid); + + copula.Theta = -0.5; + Assert.IsTrue(copula.ParametersValid); + } + + /// + /// Test Clone produces an independent copy with deep-copied marginals. + /// + [TestMethod] + public void Test_Clone() + { + var copula = new AMHCopula(0.4, new Normal(100, 10), new Gumbel(50, 5)); + var clone = copula.Clone() as AMHCopula; + Assert.IsNotNull(clone); + Assert.AreEqual(copula.Theta, clone.Theta); + Assert.IsTrue(clone.ParametersValid); + + // Marginals are deep-copied (distributions memoize lazily, so clones must + // not share marginal instances), with identical quantiles + Assert.AreNotSame(copula.MarginalDistributionX, clone.MarginalDistributionX); + Assert.AreNotSame(copula.MarginalDistributionY, clone.MarginalDistributionY); + Assert.AreEqual(copula.MarginalDistributionY.InverseCDF(0.9), clone.MarginalDistributionY.InverseCDF(0.9)); + + // Mutating the clone should not affect the original + clone.Theta = 0.8; + Assert.AreEqual(0.4, copula.Theta); + } } } diff --git a/Test_Numerics/Distributions/Bivariate Copulas/Test_ClaytonCopula.cs b/Test_Numerics/Distributions/Bivariate Copulas/Test_ClaytonCopula.cs index b112dc5a..add4d080 100644 --- a/Test_Numerics/Distributions/Bivariate Copulas/Test_ClaytonCopula.cs +++ b/Test_Numerics/Distributions/Bivariate Copulas/Test_ClaytonCopula.cs @@ -183,5 +183,46 @@ public void Test_TailDependence() var copulaZero = new ClaytonCopula(0.0); Assert.AreEqual(0.0, copulaZero.LowerTailDependence, 1E-10, "λ_L should be 0 for θ ≤ 0."); } + + /// + /// Test that ParametersValid tracks the dependency parameter's valid range. + /// Regression: the Archimedean base ValidateParameter returned a non-null + /// sentinel for valid parameters, leaving ParametersValid permanently false. + /// + [TestMethod] + public void Test_ParametersValid() + { + var copula = new ClaytonCopula(2.0); + Assert.IsTrue(copula.ParametersValid); + + copula.Theta = -2.0; // below the minimum of -1 + Assert.IsFalse(copula.ParametersValid); + + copula.Theta = 3.0; + Assert.IsTrue(copula.ParametersValid); + } + + /// + /// Test Clone produces an independent copy with deep-copied marginals. + /// + [TestMethod] + public void Test_Clone() + { + var copula = new ClaytonCopula(2.0, new Normal(100, 10), new Gumbel(50, 5)); + var clone = copula.Clone() as ClaytonCopula; + Assert.IsNotNull(clone); + Assert.AreEqual(copula.Theta, clone.Theta); + Assert.IsTrue(clone.ParametersValid); + + // Marginals are deep-copied (distributions memoize lazily, so clones must + // not share marginal instances), with identical quantiles + Assert.AreNotSame(copula.MarginalDistributionX, clone.MarginalDistributionX); + Assert.AreNotSame(copula.MarginalDistributionY, clone.MarginalDistributionY); + Assert.AreEqual(copula.MarginalDistributionY.InverseCDF(0.9), clone.MarginalDistributionY.InverseCDF(0.9)); + + // Mutating the clone should not affect the original + clone.Theta = 4.0; + Assert.AreEqual(2.0, copula.Theta); + } } } diff --git a/Test_Numerics/Distributions/Bivariate Copulas/Test_FrankCopula.cs b/Test_Numerics/Distributions/Bivariate Copulas/Test_FrankCopula.cs index 81348000..0362eec2 100644 --- a/Test_Numerics/Distributions/Bivariate Copulas/Test_FrankCopula.cs +++ b/Test_Numerics/Distributions/Bivariate Copulas/Test_FrankCopula.cs @@ -161,5 +161,38 @@ public void Test_TailDependence() Assert.AreEqual(0.0, copulaNeg.UpperTailDependence, 1E-10, "Frank copula should have no upper tail dependence for negative θ."); Assert.AreEqual(0.0, copulaNeg.LowerTailDependence, 1E-10, "Frank copula should have no lower tail dependence for negative θ."); } + + /// + /// Test that ParametersValid holds across the Frank copula's unbounded θ range. + /// + [TestMethod] + public void Test_ParametersValid() + { + Assert.IsTrue(new FrankCopula(4.2).ParametersValid); + Assert.IsTrue(new FrankCopula(-7.5).ParametersValid); + } + + /// + /// Test Clone produces an independent copy with deep-copied marginals. + /// + [TestMethod] + public void Test_Clone() + { + var copula = new FrankCopula(4.2, new Normal(100, 10), new Gumbel(50, 5)); + var clone = copula.Clone() as FrankCopula; + Assert.IsNotNull(clone); + Assert.AreEqual(copula.Theta, clone.Theta); + Assert.IsTrue(clone.ParametersValid); + + // Marginals are deep-copied (distributions memoize lazily, so clones must + // not share marginal instances), with identical quantiles + Assert.AreNotSame(copula.MarginalDistributionX, clone.MarginalDistributionX); + Assert.AreNotSame(copula.MarginalDistributionY, clone.MarginalDistributionY); + Assert.AreEqual(copula.MarginalDistributionY.InverseCDF(0.9), clone.MarginalDistributionY.InverseCDF(0.9)); + + // Mutating the clone should not affect the original + clone.Theta = 2.0; + Assert.AreEqual(4.2, copula.Theta); + } } } diff --git a/Test_Numerics/Distributions/Bivariate Copulas/Test_GumbelCopula.cs b/Test_Numerics/Distributions/Bivariate Copulas/Test_GumbelCopula.cs index ac96ed4f..210e5c58 100644 --- a/Test_Numerics/Distributions/Bivariate Copulas/Test_GumbelCopula.cs +++ b/Test_Numerics/Distributions/Bivariate Copulas/Test_GumbelCopula.cs @@ -184,5 +184,46 @@ public void Test_TailDependence() Assert.IsGreaterThan(0.99, copulaLarge.UpperTailDependence, "λ_U should approach 1 for large θ."); } + /// + /// Test that ParametersValid tracks the dependency parameter's valid range. + /// Regression: the Archimedean base ValidateParameter returned a non-null + /// sentinel for valid parameters, leaving ParametersValid permanently false. + /// + [TestMethod] + public void Test_ParametersValid() + { + var copula = new GumbelCopula(2.0); + Assert.IsTrue(copula.ParametersValid); + + copula.Theta = 0.5; // below the minimum of 1 + Assert.IsFalse(copula.ParametersValid); + + copula.Theta = 3.0; + Assert.IsTrue(copula.ParametersValid); + } + + /// + /// Test Clone produces an independent copy with deep-copied marginals. + /// + [TestMethod] + public void Test_Clone() + { + var copula = new GumbelCopula(2.0, new Normal(100, 10), new Gumbel(50, 5)); + var clone = copula.Clone() as GumbelCopula; + Assert.IsNotNull(clone); + Assert.AreEqual(copula.Theta, clone.Theta); + Assert.IsTrue(clone.ParametersValid); + + // Marginals are deep-copied (distributions memoize lazily, so clones must + // not share marginal instances), with identical quantiles + Assert.AreNotSame(copula.MarginalDistributionX, clone.MarginalDistributionX); + Assert.AreNotSame(copula.MarginalDistributionY, clone.MarginalDistributionY); + Assert.AreEqual(copula.MarginalDistributionY.InverseCDF(0.9), clone.MarginalDistributionY.InverseCDF(0.9)); + + // Mutating the clone should not affect the original + clone.Theta = 3.0; + Assert.AreEqual(2.0, copula.Theta); + } + } } diff --git a/Test_Numerics/Distributions/Bivariate Copulas/Test_JoeCopula.cs b/Test_Numerics/Distributions/Bivariate Copulas/Test_JoeCopula.cs index 2fd3aeb8..d2dc3e19 100644 --- a/Test_Numerics/Distributions/Bivariate Copulas/Test_JoeCopula.cs +++ b/Test_Numerics/Distributions/Bivariate Copulas/Test_JoeCopula.cs @@ -173,5 +173,46 @@ public void Test_TailDependence() Assert.IsGreaterThan(0.99, copulaLarge.UpperTailDependence, "λ_U should approach 1 for large θ."); } + /// + /// Test that ParametersValid tracks the dependency parameter's valid range. + /// Regression: the Archimedean base ValidateParameter returned a non-null + /// sentinel for valid parameters, leaving ParametersValid permanently false. + /// + [TestMethod] + public void Test_ParametersValid() + { + var copula = new JoeCopula(2.0); + Assert.IsTrue(copula.ParametersValid); + + copula.Theta = 0.5; // below the minimum of 1 + Assert.IsFalse(copula.ParametersValid); + + copula.Theta = 3.0; + Assert.IsTrue(copula.ParametersValid); + } + + /// + /// Test Clone produces an independent copy with deep-copied marginals. + /// + [TestMethod] + public void Test_Clone() + { + var copula = new JoeCopula(2.0, new Normal(100, 10), new Gumbel(50, 5)); + var clone = copula.Clone() as JoeCopula; + Assert.IsNotNull(clone); + Assert.AreEqual(copula.Theta, clone.Theta); + Assert.IsTrue(clone.ParametersValid); + + // Marginals are deep-copied (distributions memoize lazily, so clones must + // not share marginal instances), with identical quantiles + Assert.AreNotSame(copula.MarginalDistributionX, clone.MarginalDistributionX); + Assert.AreNotSame(copula.MarginalDistributionY, clone.MarginalDistributionY); + Assert.AreEqual(copula.MarginalDistributionY.InverseCDF(0.9), clone.MarginalDistributionY.InverseCDF(0.9)); + + // Mutating the clone should not affect the original + clone.Theta = 3.0; + Assert.AreEqual(2.0, copula.Theta); + } + } } diff --git a/Test_Numerics/Distributions/Bivariate Copulas/Test_NormalCopula.cs b/Test_Numerics/Distributions/Bivariate Copulas/Test_NormalCopula.cs index 6e43b5c0..56ac6357 100644 --- a/Test_Numerics/Distributions/Bivariate Copulas/Test_NormalCopula.cs +++ b/Test_Numerics/Distributions/Bivariate Copulas/Test_NormalCopula.cs @@ -205,5 +205,44 @@ public void Test_TailDependence() Assert.AreEqual(0.0, copulaHigh.UpperTailDependence, 1E-10, "Normal copula should have no tail dependence even with high ρ."); Assert.AreEqual(0.0, copulaHigh.LowerTailDependence, 1E-10, "Normal copula should have no tail dependence even with high ρ."); } + + /// + /// Test that ParametersValid tracks the correlation parameter's valid range. + /// + [TestMethod] + public void Test_ParametersValid() + { + var copula = new NormalCopula(0.5); + Assert.IsTrue(copula.ParametersValid); + + copula.Theta = 1.5; // above the maximum of 1 + Assert.IsFalse(copula.ParametersValid); + + copula.Theta = -0.35; + Assert.IsTrue(copula.ParametersValid); + } + + /// + /// Test Clone produces an independent copy with deep-copied marginals. + /// + [TestMethod] + public void Test_Clone() + { + var copula = new NormalCopula(0.5, new Normal(100, 10), new Gumbel(50, 5)); + var clone = copula.Clone() as NormalCopula; + Assert.IsNotNull(clone); + Assert.AreEqual(copula.Theta, clone.Theta); + Assert.IsTrue(clone.ParametersValid); + + // Marginals are deep-copied (distributions memoize lazily, so clones must + // not share marginal instances), with identical quantiles + Assert.AreNotSame(copula.MarginalDistributionX, clone.MarginalDistributionX); + Assert.AreNotSame(copula.MarginalDistributionY, clone.MarginalDistributionY); + Assert.AreEqual(copula.MarginalDistributionY.InverseCDF(0.9), clone.MarginalDistributionY.InverseCDF(0.9)); + + // Mutating the clone should not affect the original + clone.Theta = 0.9; + Assert.AreEqual(0.5, copula.Theta); + } } } diff --git a/Test_Numerics/Distributions/Bivariate Copulas/Test_StudentTCopula.cs b/Test_Numerics/Distributions/Bivariate Copulas/Test_StudentTCopula.cs index 7ae70a2e..4b91d9a8 100644 --- a/Test_Numerics/Distributions/Bivariate Copulas/Test_StudentTCopula.cs +++ b/Test_Numerics/Distributions/Bivariate Copulas/Test_StudentTCopula.cs @@ -353,6 +353,15 @@ public void Test_Clone() // Mutating clone should not affect original clone.Theta = 0.1; Assert.AreEqual(0.5, copula.Theta); + + // Marginals are deep-copied (distributions memoize lazily, so clones must + // not share marginal instances), with identical quantiles + var withMarginals = new StudentTCopula(0.5, 10, new Normal(100, 10), new Gumbel(50, 5)); + var marginalClone = withMarginals.Clone() as StudentTCopula; + Assert.IsNotNull(marginalClone); + Assert.AreNotSame(withMarginals.MarginalDistributionX, marginalClone.MarginalDistributionX); + Assert.AreNotSame(withMarginals.MarginalDistributionY, marginalClone.MarginalDistributionY); + Assert.AreEqual(withMarginals.MarginalDistributionY.InverseCDF(0.9), marginalClone.MarginalDistributionY.InverseCDF(0.9)); } private double[] data1 = new double[] { 122.094066003419, 92.8321267206161, 86.4920318705377, 87.6183663113541, 102.558777787492, 103.627475117762, 127.084948716539, 105.908684131013, 110.065795957654, 105.924647125867, 110.009738155469, 126.490833800772, 64.1264871206211, 81.3150800229481, 92.0780134395721, 106.040322550555, 113.158086143066, 117.051057784044, 127.110531266645, 108.907371862136, 105.476247114194, 108.629403495407, 98.7803988364997, 93.217925588845, 97.7219451830075, 109.178093756809, 137.69504856252, 106.884615327674, 112.139177456202, 85.7416217661797, 71.0610938629716, 112.644166631765, 119.545871678548, 70.5169833274982, 99.6896817997206, 100.987892854545, 103.659280253554, 75.6075621013066, 118.810868919796, 109.113664695226, 113.636425353944, 100.008375355612, 113.178917359795, 80.4269472604342, 88.3638384448237, 90.2905074656314, 98.7995143316863, 98.4698060067802, 108.279297570816, 86.1578437055905, 101.183725242941, 85.5531148952956, 111.024195253862, 121.934506174556, 104.169993666179, 84.4652994609478, 99.6099259747033, 95.3130792386208, 115.45680252817, 120.213139478586, 95.5691788140058, 92.7950300448044, 102.58430893827, 86.7105161576407, 82.8059368562185, 107.335705516294, 112.603259240932, 102.780778760832, 128.958090528336, 105.139162595628, 118.272661482198, 99.8275937885748, 94.2856024560543, 108.48679008009, 100.147734981682, 88.7006383425785, 89.6441478272035, 112.24266306884, 99.8184811468069, 120.592090049738, 124.023170133661, 101.250961381805, 90.0000027551006, 108.781064635426, 94.9203320035987, 99.9491821782837, 88.7473944659517, 94.3643253649856, 105.814317118952, 92.6866900633813, 111.020330544613, 111.676189456988, 115.70235103978, 124.659106152655, 81.3866270495082, 120.178528245778, 93.6511977805724, 114.099368762143, 119.062045395294, 74.1998497412903 }; diff --git a/Test_Numerics/Distributions/Multivariate/Test_BivariateEmpirical.cs b/Test_Numerics/Distributions/Multivariate/Test_BivariateEmpirical.cs index c42343e8..da0646ba 100644 --- a/Test_Numerics/Distributions/Multivariate/Test_BivariateEmpirical.cs +++ b/Test_Numerics/Distributions/Multivariate/Test_BivariateEmpirical.cs @@ -1,4 +1,5 @@ using Microsoft.VisualStudio.TestTools.UnitTesting; +using System; using Numerics.Distributions; namespace Distributions.Multivariate @@ -49,8 +50,31 @@ public void Test_BivariateEmp() double true_cdf5 = 0.386806419d; bv.ProbabilityTransform = Numerics.Data.Transform.NormalZ; double cdf5 = bv.CDF(3.22d, 8.15d); + Assert.AreEqual(cdf5, true_cdf5, 1E-6); } + /// + /// Verifies that replacing the parameter grid invalidates the cached transformed interpolator. + /// + [TestMethod] + public void SetParametersInvalidatesBilinearInterpolator() + { + var x1 = new[] { 1d, 10d }; + var x2 = new[] { 1d, 10d }; + var original = new[,] { { 0.1d, 0.2d }, { 0.3d, 0.4d } }; + var replacement = new[,] { { 0.6d, 0.7d }, { 0.8d, 0.9d } }; + var distribution = new BivariateEmpirical(x1, x2, original) + { + X1Transform = Numerics.Data.Transform.Logarithmic, + X2Transform = Numerics.Data.Transform.Logarithmic + }; + + double coordinate = Math.Sqrt(10d); + Assert.AreEqual(0.25d, distribution.CDF(coordinate, coordinate), 1E-12); + distribution.SetParameters(x1, x2, replacement); + Assert.AreEqual(0.75d, distribution.CDF(coordinate, coordinate), 1E-12); + } + } } diff --git a/Test_Numerics/Distributions/Multivariate/Test_Dirichlet.cs b/Test_Numerics/Distributions/Multivariate/Test_Dirichlet.cs index b8f5321d..353609d0 100644 --- a/Test_Numerics/Distributions/Multivariate/Test_Dirichlet.cs +++ b/Test_Numerics/Distributions/Multivariate/Test_Dirichlet.cs @@ -132,7 +132,7 @@ public void Test_Mode() } /// - /// Test that Mode throws when any alpha <= 1. + /// Test that Mode throws when any alpha <= 1. /// [TestMethod] public void Test_Mode_Invalid() diff --git a/Test_Numerics/Distributions/Multivariate/Test_MultivariateNormal.cs b/Test_Numerics/Distributions/Multivariate/Test_MultivariateNormal.cs index 93d227a5..b67b5dd8 100644 --- a/Test_Numerics/Distributions/Multivariate/Test_MultivariateNormal.cs +++ b/Test_Numerics/Distributions/Multivariate/Test_MultivariateNormal.cs @@ -20,6 +20,129 @@ namespace Distributions.Multivariate public class Test_MultivariateNormal { + /// + /// Asserts the action throws an ArgumentOutOfRangeException (net481-compatible). + /// + /// The action expected to throw. + private static void AssertThrowsOutOfRange(Action action) + { + try + { + action(); + } + catch (ArgumentOutOfRangeException) + { + return; + } + Assert.Fail("Expected an ArgumentOutOfRangeException."); + } + + + /// + /// Verifies the non-throwing mutable-covariance path: a positive-definite swap + /// evaluates normally, a non-positive-definite swap marks the density invalid + /// (LogPDF = -inf, PDF = 0) without throwing, and a subsequent valid swap + /// restores evaluation. + /// + [TestMethod()] + public void Test_TrySetCovariance_NonThrowingInvalidState() + { + var mvn = new MultivariateNormal(new[] { 1d, 2d }, new[,] { { 1d, 0.3d }, { 0.3d, 2d } }); + double validLogPdf = mvn.LogPDF(new[] { 1.2d, 1.8d }); + Assert.IsTrue(mvn.IsDensityValid); + Assert.IsTrue(!double.IsNaN(validLogPdf) && !double.IsInfinity(validLogPdf)); + + // A non-positive-definite covariance: correlation beyond one. + Assert.IsFalse(mvn.TrySetCovariance(new[,] { { 1d, 1.5d }, { 1.5d, 1d } })); + Assert.IsFalse(mvn.IsDensityValid); + Assert.IsTrue(double.IsNegativeInfinity(mvn.LogPDF(new[] { 1.2d, 1.8d }))); + Assert.AreEqual(0d, mvn.PDF(new[] { 1.2d, 1.8d }), 0d); + + // Restore with the original covariance: the density comes back exactly. + Assert.IsTrue(mvn.TrySetCovariance(new[,] { { 1d, 0.3d }, { 0.3d, 2d } })); + Assert.IsTrue(mvn.IsDensityValid); + Assert.AreEqual(validLogPdf, mvn.LogPDF(new[] { 1.2d, 1.8d }), 1E-12); + + // TrySetParameters also moves the mean. + Assert.IsTrue(mvn.TrySetParameters(new[] { 0d, 0d }, new[,] { { 1d, 0d }, { 0d, 1d } })); + Assert.AreEqual(-Math.Log(2d * Math.PI), mvn.LogPDF(new[] { 0d, 0d }), 1E-12); + } + + /// + /// Verifies the marginal utility: the sub-mean and sub-covariance at the kept + /// indices, including a reordered subset, and the index validation matrix. + /// + [TestMethod()] + public void Test_Marginal_SubsetAndValidation() + { + var mean = new[] { 1d, 2d, 3d }; + var covariance = new[,] + { + { 4.0d, 1.2d, 0.5d }, + { 1.2d, 9.0d, 2.1d }, + { 0.5d, 2.1d, 16.0d } + }; + var mvn = new MultivariateNormal(mean, covariance); + + var marginal = mvn.Marginal(2, 0); + Assert.AreEqual(2, marginal.Dimension); + Assert.AreEqual(3d, marginal.Mean[0], 0d); + Assert.AreEqual(1d, marginal.Mean[1], 0d); + Assert.AreEqual(16d, marginal.Covariance[0, 0], 0d); + Assert.AreEqual(0.5d, marginal.Covariance[0, 1], 0d); + Assert.AreEqual(4d, marginal.Covariance[1, 1], 0d); + + // The marginal of a joint Gaussian is its sub-Gaussian: densities agree + // with a directly constructed distribution. + var direct = new MultivariateNormal(new[] { 3d, 1d }, + new[,] { { 16d, 0.5d }, { 0.5d, 4d } }); + Assert.AreEqual(direct.LogPDF(new[] { 2.5d, 1.4d }), marginal.LogPDF(new[] { 2.5d, 1.4d }), 1E-12); + + AssertThrowsOutOfRange(() => mvn.Marginal()); + AssertThrowsOutOfRange(() => mvn.Marginal(0, 0)); + AssertThrowsOutOfRange(() => mvn.Marginal(3)); + } + + /// + /// Verifies the conditional utility against hand-computed Gaussian + /// conditioning: for the bivariate case, X1 | X2 = x2 has mean + /// mu1 + rho*(s1/s2)*(x2 - mu2) and variance s1^2*(1 - rho^2); a trivariate + /// case checks the Schur complement, and observing everything throws. + /// + [TestMethod()] + public void Test_Conditional_ClosedFormsAndValidation() + { + // Bivariate closed form. + double mu1 = 4d, mu2 = 2d, s1 = 2d, s2 = 3d, rho = 0.6d; + var bivariate = MultivariateNormal.Bivariate(mu1, mu2, s1, s2, rho); + var conditional = bivariate.Conditional(new[] { 1 }, new[] { 3.5d }); + Assert.AreEqual(1, conditional.Dimension); + Assert.AreEqual(mu1 + rho * (s1 / s2) * (3.5d - mu2), conditional.Mean[0], 1E-12); + Assert.AreEqual(s1 * s1 * (1d - rho * rho), conditional.Covariance[0, 0], 1E-12); + + // Trivariate Schur complement, observing the middle dimension. + var mean = new[] { 1d, 2d, 3d }; + var covariance = new[,] + { + { 4.0d, 1.2d, 0.5d }, + { 1.2d, 9.0d, 2.1d }, + { 0.5d, 2.1d, 16.0d } + }; + var trivariate = new MultivariateNormal(mean, covariance); + var given = trivariate.Conditional(new[] { 1 }, new[] { 4.0d }); + // mu_c = mu_{0,2} + Sigma_{[0,2],1} * (4 - 2) / 9 + Assert.AreEqual(1d + 1.2d * 2d / 9d, given.Mean[0], 1E-12); + Assert.AreEqual(3d + 2.1d * 2d / 9d, given.Mean[1], 1E-12); + // Sigma_c = Sigma_{[0,2],[0,2]} - outer(Sigma_{[0,2],1}) / 9 + Assert.AreEqual(4d - 1.2d * 1.2d / 9d, given.Covariance[0, 0], 1E-12); + Assert.AreEqual(0.5d - 1.2d * 2.1d / 9d, given.Covariance[0, 1], 1E-12); + Assert.AreEqual(16d - 2.1d * 2.1d / 9d, given.Covariance[1, 1], 1E-12); + + AssertThrowsOutOfRange(() => trivariate.Conditional(new[] { 0, 1, 2 }, new[] { 1d, 2d, 3d })); + AssertThrowsOutOfRange(() => trivariate.Conditional(new[] { 1 }, new[] { 1d, 2d })); + } + + /// /// Verified using Accord.Net /// @@ -108,6 +231,142 @@ public void Test_MultivariateNormalCDF_Fortran() } + /// + /// Verifies the documented invalid-dimension termination status without entering MVNDNT initialization. + /// + [TestMethod] + public void Test_MVNDST_InvalidDimensionReturnsStatusTwo() + { + var distribution = new MultivariateNormal(1); + double error = 0d; + double value = 0d; + int inform = 0; + + distribution.MVNDST(0, new double[0], new double[0], new int[0], new double[0], 1, 0d, 0d, ref error, ref value, ref inform); + + Assert.AreEqual(2, inform); + Assert.AreEqual(0d, value, 0d); + Assert.AreEqual(1d, error, 0d); + } + + /// + /// Verifies that dimensions flagged as completely unbounded collapse to probability one with successful status. + /// + [TestMethod] + public void Test_MVNDST_AllUnboundedReturnsExactProbability() + { + var distribution = new MultivariateNormal(2); + double error = -1d; + double value = -1d; + int inform = -1; + + distribution.MVNDST( + 2, + new[] { 0d, 0d }, + new[] { 0d, 0d }, + new[] { -1, -1 }, + new[] { 0d }, + 1, + 0d, + 0d, + ref error, + ref value, + ref inform); + + Assert.AreEqual(0, inform); + Assert.AreEqual(1d, value, 0d); + Assert.AreEqual(0d, error, 0d); + } + + /// + /// Verifies the exact one-active-dimension Normal limit and its successful initialization status. + /// + [TestMethod] + public void Test_MVNDST_OneActiveDimensionMatchesNormalProbability() + { + var distribution = new MultivariateNormal(2); + double error = 0d; + double value = 0d; + int inform = -1; + + distribution.MVNDST( + 2, + new[] { 0d, 0d }, + new[] { 0d, 0d }, + new[] { 0, -1 }, + new[] { 0.7d }, + 1, + 0d, + 0d, + ref error, + ref value, + ref inform); + + Assert.AreEqual(0, inform); + Assert.AreEqual(0.5d, value, 1E-15); + Assert.AreEqual(2E-16, error, 0d); + } + + /// + /// Verifies the analytical bivariate collapse against the zero-threshold quadrant probability for correlation one-half. + /// + [TestMethod] + public void Test_MVNDST_BivariateCollapseMatchesAnalyticalProbability() + { + var distribution = new MultivariateNormal(2); + double error = 0d; + double value = 0d; + int inform = -1; + + distribution.MVNDST( + 2, + new[] { 0d, 0d }, + new[] { 0d, 0d }, + new[] { 0, 0 }, + new[] { 0.5d }, + 1, + 0d, + 0d, + ref error, + ref value, + ref inform); + + Assert.AreEqual(0, inform); + Assert.AreEqual(1d / 3d, value, 1E-12); + Assert.AreEqual(2E-16, error, 0d); + } + + /// + /// Verifies that only the lattice integration stage changes successful initialization into a budget-exhaustion status. + /// + [TestMethod] + public void Test_MVNDST_InsufficientBudgetReturnsStatusOne() + { + var distribution = new MultivariateNormal(3) { MVNUNI = new MersenneTwister(12345) }; + double error = 0d; + double value = 0d; + int inform = -1; + + distribution.MVNDST( + 3, + new[] { 0d, 0d, 0d }, + new[] { 0.2d, 0.5d, 1d }, + new[] { 0, 0, 0 }, + new[] { 0.25d, 0.1d, 0.3d }, + 1, + 0d, + 0d, + ref error, + ref value, + ref inform); + + Assert.AreEqual(1, inform); + Assert.IsTrue(value >= 0d && value <= 1d); + Assert.IsGreaterThanOrEqualTo(0d, error); + Assert.IsFalse(double.IsNaN(error)); + Assert.IsFalse(double.IsInfinity(error)); + } + /// /// Verified against R "mvtnorm" package /// @@ -231,5 +490,92 @@ public void Test_MultivariateNormalCDF_R_Independent() } - } + /// + /// Verifies that perfectly correlated variables collapse to their tightest univariate bound. + /// + [TestMethod] + public void Test_CDF_PerfectCorrelationCollapsesToNormal() + { + double probability = EvaluateStandardCdf(new[] { 0.8d, -0.3d, 0.2d }, new[] { 1d, 1d, 1d }); + Assert.AreEqual(Normal.StandardCDF(-0.3d), probability, 1E-10); + } + + /// + /// Verifies that perfectly anticorrelated variables collapse to a bounded normal interval. + /// + [TestMethod] + public void Test_CDF_PerfectAnticorrelationCollapsesToNormalInterval() + { + double probability = EvaluateStandardCdf(new[] { 0.7d, 0.2d }, new[] { -1d }); + double expected = Normal.StandardCDF(0.7d) - Normal.StandardCDF(-0.2d); + + Assert.AreEqual(expected, probability, 1E-10); + } + + /// + /// Verifies equivalent results for permutations of a rank-deficient covariance matrix. + /// + [TestMethod] + public void Test_CDF_PermutedRankDeficientMatricesCollapseAnalytically() + { + double expected = Normal.StandardCDF(0.2d) * Normal.StandardCDF(-0.4d); + Assert.AreEqual( + expected, + EvaluateStandardCdf(new[] { 0.2d, -0.4d, 0.7d }, new[] { 0d, 1d, 0d }), + 1E-10); + Assert.AreEqual( + expected, + EvaluateStandardCdf(new[] { 0.7d, 0.2d, -0.4d }, new[] { 1d, 0d, 0d }), + 1E-10); + } + + /// + /// Verifies the unchanged nonsingular path immediately above perfect correlation. + /// + [TestMethod] + public void Test_CDF_NearSingularCorrelationMatchesBivariateFormula() + { + const double correlation = 1d - 1E-12; + var distribution = MultivariateNormal.Bivariate(0d, 0d, 1d, 1d, correlation); + distribution.MVNUNI = new MersenneTwister(12345); + double expected = 0.25d + Math.Asin(correlation) / (2d * Math.PI); + + Assert.AreEqual(expected, distribution.CDF(new[] { 0d, 0d }), 1E-10); + } + + /// + /// Evaluates a standard multivariate-normal CDF through the public Genz integration entry point. + /// + /// The upper integration limits. + /// The packed strict-lower-triangle correlation coefficients. + /// The evaluated multivariate-normal probability. + private static double EvaluateStandardCdf(double[] upper, double[] correlations) + { + int dimensions = upper.Length; + var distribution = new MultivariateNormal(dimensions) + { + MVNUNI = new MersenneTwister(12345) + }; + var lower = new double[dimensions]; + var infinities = new int[dimensions]; + double error = 0d; + double value = 0d; + int inform = 0; + + distribution.MVNDST( + dimensions, + lower, + upper, + infinities, + correlations, + 25000, + 1E-10, + 0d, + ref error, + ref value, + ref inform); + return value; + } + } + } diff --git a/Test_Numerics/Distributions/Multivariate/Test_MultivariateStudentT.cs b/Test_Numerics/Distributions/Multivariate/Test_MultivariateStudentT.cs index 7210e8dc..3a96a283 100644 --- a/Test_Numerics/Distributions/Multivariate/Test_MultivariateStudentT.cs +++ b/Test_Numerics/Distributions/Multivariate/Test_MultivariateStudentT.cs @@ -204,7 +204,7 @@ public void Test_PDF_ReducesToUnivariate() { double mvtPdf = mvt.PDF(new[] { x }); double univPdf = univT.PDF(x); - Assert.AreEqual(univPdf / sigma, mvtPdf, 1E-12, + Assert.AreEqual(univPdf, mvtPdf, 1E-12, $"1D MVT PDF does not match univariate StudentT PDF at x={x}"); } } diff --git a/Test_Numerics/Distributions/Test_ParameterValidity.cs b/Test_Numerics/Distributions/Test_ParameterValidity.cs new file mode 100644 index 00000000..46fa3ab0 --- /dev/null +++ b/Test_Numerics/Distributions/Test_ParameterValidity.cs @@ -0,0 +1,225 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using Microsoft.VisualStudio.TestTools.UnitTesting; +using Numerics.Distributions; +using Numerics.Distributions.Copulas; + +namespace Distributions +{ + /// + /// Regression tests for parameter-validity state across univariate, multivariate, and bivariate distributions. + /// + [TestClass] + public class Test_ParameterValidity + { + /// + /// Verifies every factory-supported flat univariate parameter vector transitions from valid to invalid and back. + /// + [TestMethod] + public void UnivariateFlatSettersTrackFinalParameterTuple() + { + foreach (UnivariateDistributionType type in Enum.GetValues(typeof(UnivariateDistributionType))) + { + if (!SupportsFlatParameterReplacement(type)) continue; + + var distribution = UnivariateDistributionFactory.CreateDistribution(type); + Assert.IsTrue(distribution.ParametersValid, $"{type} default parameters must be valid."); + double[] validParameters = distribution.GetParameters.ToArray(); + double[] invalidFirst = validParameters.ToArray(); + invalidFirst[0] = double.NaN; + distribution.SetParameters(invalidFirst); + Assert.IsFalse(distribution.ParametersValid, $"{type} must reject NaN in its first parameter."); + + invalidFirst[0] = double.PositiveInfinity; + distribution.SetParameters(invalidFirst); + Assert.IsFalse(distribution.ParametersValid, $"{type} must reject infinity in its first parameter."); + distribution.SetParameters(validParameters); + Assert.IsTrue(distribution.ParametersValid, $"{type} must recover after restoring its first parameter."); + + double[] invalidFinal = validParameters.ToArray(); + invalidFinal[invalidFinal.Length - 1] = double.NaN; + distribution.SetParameters(invalidFinal); + Assert.IsFalse(distribution.ParametersValid, $"{type} must reject NaN in its final parameter."); + + invalidFinal[invalidFinal.Length - 1] = double.PositiveInfinity; + distribution.SetParameters(invalidFinal); + Assert.IsFalse(distribution.ParametersValid, $"{type} must reject infinity in its final parameter."); + distribution.SetParameters(validParameters); + Assert.IsTrue(distribution.ParametersValid, $"{type} must recover after restoring its final parameter."); + CollectionAssert.AreEqual(validParameters, distribution.GetParameters, $"{type} must store the restored tuple."); + } + } + + /// + /// Verifies data-backed and composite univariate distributions refresh validity after every replacement path. + /// + [TestMethod] + public void SpecializedUnivariateSettersTrackFinalState() + { + var empirical = new EmpiricalDistribution(new[] { 0d, 1d }, new[] { 0d, 1d }); + empirical.SetParameters(new[] { 0d, 1d }, new[] { 0d, 2d }); + Assert.IsFalse(empirical.ParametersValid); + Assert.Throws(() => empirical.CDF(0.5d)); + Assert.Throws(() => empirical.SetParameters(new[] { 0d }, new[] { 0d, 1d })); + empirical.SetParameters(new[] { 0d, 1d }, new[] { 0d, 1d }); + Assert.IsTrue(empirical.ParametersValid); + empirical.SetParameters(new[] { 0d, 0d, 1d }, new[] { 0d, 0.5d, 1d }); + Assert.IsTrue(empirical.ParametersValid); + empirical.SetParameters(new[] { 0d, -1d }, new[] { 0d, 1d }); + Assert.IsFalse(empirical.ParametersValid); + empirical.SetParameters(new[] { 0d, 1d }, new[] { 0d, 1d }); + Assert.IsTrue(empirical.ParametersValid); + + + var kernel = new KernelDensity(new[] { -1d, 0d, 1d }); + kernel.Bandwidth = double.NaN; + Assert.IsFalse(kernel.ParametersValid); + kernel.Bandwidth = double.PositiveInfinity; + Assert.IsFalse(kernel.ParametersValid); + kernel.Bandwidth = 0.5d; + Assert.IsTrue(kernel.ParametersValid); + + var mixture = new Mixture( + new[] { 0.8d, 0.8d }, + new UnivariateDistributionBase[] { new Normal(), new Exponential() }); + Assert.IsFalse(mixture.ParametersValid); + mixture.SetParameters( + new[] { 0.5d, 0.5d }, + new UnivariateDistributionBase[] { new Normal(), new Exponential() }); + Assert.IsTrue(mixture.ParametersValid); + mixture.SetParameters(new[] { 0.5d, 0.5d }, new[] { 0d, -1d, 0d, 1d }); + Assert.IsFalse(mixture.ParametersValid); + + var emptyMixture = new Mixture(Array.Empty(), Array.Empty()); + Assert.IsFalse(emptyMixture.ParametersValid); + + var competingRisks = new CompetingRisks(Array.Empty()); + var zeroInflatedMixture = new Mixture( + new[] { 0.4d, 0.5d }, + new UnivariateDistributionBase[] { new Normal(), new Exponential() }); + zeroInflatedMixture.IsZeroInflated = true; + zeroInflatedMixture.ZeroWeight = 0.1d; + Assert.IsTrue(zeroInflatedMixture.ParametersValid); + zeroInflatedMixture.IsZeroInflated = false; + Assert.IsFalse(zeroInflatedMixture.ParametersValid); + Assert.IsFalse(competingRisks.ParametersValid); + competingRisks.SetParameters(new UnivariateDistributionBase[] { new Normal(0d, -1d) }); + Assert.IsFalse(competingRisks.ParametersValid); + competingRisks.SetParameters(new UnivariateDistributionBase[] { new Normal() }); + Assert.IsTrue(competingRisks.ParametersValid); + Assert.Throws(() => competingRisks.SetParameters(new[] { 0d })); + + var truncated = new TruncatedDistribution(new Normal(), -1d, 1d); + Assert.IsTrue(truncated.ParametersValid); + truncated.SetParameters(new[] { 0d, 1d, 2d, 1d }); + Assert.IsFalse(truncated.ParametersValid); + Assert.Throws(() => truncated.CDF(1.5d)); + truncated.SetParameters(new[] { 2d, 3d, -1d, 1d }); + Assert.IsTrue(truncated.ParametersValid); + CollectionAssert.AreEqual(new[] { 2d, 3d, -1d, 1d }, truncated.GetParameters); + truncated.SetParameters(new[] { 2d, 3d, 100d, 101d }); + Assert.Throws(() => truncated.SetParameters(new[] { 0d, 1d, -1d })); + Assert.IsFalse(truncated.ParametersValid); + } + + /// + /// Verifies multivariate setters reject non-finite values before mutation and bivariate empirical data can recover. + /// + [TestMethod] + public void MultivariateSettersRejectNonFiniteValuesTransactionally() + { + var covariance = new[,] { { 1d, 0.25d }, { 0.25d, 2d } }; + var normal = new MultivariateNormal(new[] { 1d, 2d }, covariance); + Assert.Throws(() => normal.SetParameters(new[] { double.NaN, 2d }, covariance)); + CollectionAssert.AreEqual(new[] { 1d, 2d }, normal.Mean); + Assert.Throws(() => normal.SetParameters(new[] { 1d, 2d }, new[,] { { 1d, double.PositiveInfinity }, { 0.25d, 2d } })); + CollectionAssert.AreEqual(new[] { 1d, 2d }, normal.Mean); + Assert.AreEqual(0.25d, normal.Covariance[0, 1]); + Assert.Throws(() => normal.SetParameters(new[] { 1d }, covariance)); + CollectionAssert.AreEqual(new[] { 1d, 2d }, normal.Mean); + normal.SetParameters(new[] { 3d, 4d }, covariance); + CollectionAssert.AreEqual(new[] { 3d, 4d }, normal.Mean); + + var student = new MultivariateStudentT(5d, new[] { 1d, 2d }, covariance); + Assert.Throws(() => student.SetParameters(double.NaN, new[] { 0d, 0d }, covariance)); + Assert.AreEqual(5d, student.DegreesOfFreedom); + Assert.Throws(() => student.SetParameters(5d, new[] { 0d, double.PositiveInfinity }, covariance)); + Assert.AreEqual(5d, student.DegreesOfFreedom); + CollectionAssert.AreEqual(new[] { 1d, 2d }, student.Location); + Assert.Throws(() => student.SetParameters(5d, new[] { 0d }, covariance)); + CollectionAssert.AreEqual(new[] { 1d, 2d }, student.Location); + student.SetParameters(7d, new[] { 3d, 4d }, covariance); + Assert.AreEqual(7d, student.DegreesOfFreedom); + CollectionAssert.AreEqual(new[] { 3d, 4d }, student.Location); + + var probabilities = new[,] { { 0.1d, 0.2d }, { 0.3d, 0.8d } }; + var empirical = new BivariateEmpirical(new[] { 0d, 1d }, new[] { 0d, 1d }, probabilities); + empirical.SetParameters(new[] { 0d, double.NaN }, new[] { 0d, 1d }, probabilities); + Assert.IsFalse(empirical.ParametersValid); + Assert.Throws(() => empirical.CDF(0.5d, 0.5d)); + empirical.SetParameters(new[] { 0d, 1d }, new[] { 0d, double.NegativeInfinity }, probabilities); + Assert.IsFalse(empirical.ParametersValid); + empirical.SetParameters(new[] { 0d, 1d }, new[] { 0d, 1d }, new[,] { { 0.1d, double.PositiveInfinity }, { 0.3d, 0.8d } }); + Assert.IsFalse(empirical.ParametersValid); + empirical.SetParameters(new[] { 0d, 1d }, new[] { 0d, 1d }, new[,] { { 0.1d }, { 0.3d } }); + Assert.IsFalse(empirical.ParametersValid); + empirical.SetParameters(new[] { 0d, 1d }, new[] { 0d, 1d }, probabilities); + Assert.IsTrue(empirical.ParametersValid); + + Assert.Throws(() => new Dirichlet(new[] { 1d, double.NaN })); + Assert.Throws(() => new Dirichlet(new[] { 1d, double.PositiveInfinity })); + Assert.Throws(() => new Multinomial(2, new[] { 0.5d, double.NaN })); + Assert.Throws(() => new Multinomial(2, new[] { 0.5d, double.PositiveInfinity })); + } + + /// + /// Verifies every copula rejects non-finite dependency parameters and recovers when restored. + /// + [TestMethod] + public void CopulasRejectNonFiniteParametersAndRecover() + { + var copulas = new BivariateCopula[] + { + new AMHCopula(0.5d), + new ClaytonCopula(2d), + new FrankCopula(4d), + new GumbelCopula(2d), + new JoeCopula(2d), + new NormalCopula(0.5d), + new StudentTCopula(0.5d, 5d) + }; + + foreach (BivariateCopula copula in copulas) + { + double validTheta = copula.Theta; + copula.Theta = double.NaN; + Assert.IsFalse(copula.ParametersValid, $"{copula.Type} must reject NaN."); + copula.Theta = double.PositiveInfinity; + Assert.IsFalse(copula.ParametersValid, $"{copula.Type} must reject infinity."); + copula.Theta = validTheta; + Assert.IsTrue(copula.ParametersValid, $"{copula.Type} must recover after restoring theta."); + } + + var student = new StudentTCopula(0.5d, 5d); + student.DegreesOfFreedom = double.PositiveInfinity; + Assert.IsFalse(student.ParametersValid); + student.DegreesOfFreedom = 5d; + Assert.IsTrue(student.ParametersValid); + } + + /// + /// Determines whether a factory-created distribution supports flattened numeric parameter replacement. + /// + /// The univariate distribution type. + /// when is supported without external components. + private static bool SupportsFlatParameterReplacement(UnivariateDistributionType type) + { + return type != UnivariateDistributionType.CompetingRisks && + type != UnivariateDistributionType.Empirical && + type != UnivariateDistributionType.KernelDensity && + type != UnivariateDistributionType.Mixture && + type != UnivariateDistributionType.UserDefined; + } + } +} diff --git a/Test_Numerics/Distributions/Univariate/Test_Beta.cs b/Test_Numerics/Distributions/Univariate/Test_Beta.cs index 24a3c40c..9d805008 100644 --- a/Test_Numerics/Distributions/Univariate/Test_Beta.cs +++ b/Test_Numerics/Distributions/Univariate/Test_Beta.cs @@ -196,6 +196,18 @@ public void Test_Mode() var b4 = new BetaDistribution(9, 1); Assert.AreEqual(1, b4.Mode); + + var lowerBoundary = new BetaDistribution(0.5d, 1.0d); + Assert.AreEqual(0.0d, lowerBoundary.Mode); + + var upperBoundary = new BetaDistribution(1.0d, 0.5d); + Assert.AreEqual(1.0d, upperBoundary.Mode); + + var uShaped = new BetaDistribution(0.5d, 0.5d); + Assert.AreEqual(0.5d, uShaped.Mode); + + var interior = new BetaDistribution(2.0d, 5.0d); + Assert.AreEqual(0.2d, interior.Mode, 1E-14); } /// diff --git a/Test_Numerics/Distributions/Univariate/Test_CompetingRisks.cs b/Test_Numerics/Distributions/Univariate/Test_CompetingRisks.cs index cbaef530..37afe4e5 100644 --- a/Test_Numerics/Distributions/Univariate/Test_CompetingRisks.cs +++ b/Test_Numerics/Distributions/Univariate/Test_CompetingRisks.cs @@ -1,5 +1,6 @@ using Microsoft.VisualStudio.TestTools.UnitTesting; using Numerics.Distributions; +using Numerics.Data.Statistics; using Numerics.Mathematics; using Numerics.Mathematics.Integration; using Numerics.Mathematics.SpecialFunctions; @@ -451,6 +452,43 @@ public void Test_LogPDF_DependentInvalidSupport_ReturnsNegativeInfinity() #region Minimum Rule - 2 Distributions + /// + /// Verifies that dependency changes rebuild the cached Gaussian copula without mutating user correlation input. + /// + [TestMethod] + public void Test_DependencyChangeInvalidatesMvnWithoutMutatingCorrelation() + { + var correlation = new[,] { { 1d, 0.5d }, { 0.5d, 1d } }; + var risks = new CompetingRisks(new IUnivariateDistribution[] + { + new Normal(), + new Normal() + }) + { + CorrelationMatrix = (double[,])correlation.Clone(), + Dependency = Probability.DependencyType.PerfectlyNegative, + MinimumOfRandomVariables = true + }; + + double perfectlyNegative = risks.CDF(0d); + + for (int row = 0; row < 2; row++) + { + for (int column = 0; column < 2; column++) + { + Assert.AreEqual(correlation[row, column], risks.CorrelationMatrix[row, column], 0d); + } + } + + risks.Dependency = Probability.DependencyType.CorrelationMatrix; + double correlated = risks.CDF(0d); + double expected = 0.75d - Math.Asin(0.5d) / (2d * Math.PI); + + Assert.AreEqual(1d, perfectlyNegative, 5E-5); + Assert.AreEqual(expected, correlated, 1E-7); + Assert.IsGreaterThan(0.1d, perfectlyNegative - correlated); + } + // Tolerances - competing risks MLE is harder than single distribution MLE private const double SHAPE_TOLERANCE_PERCENT = 0.25; // 25% relative error private const double SCALE_TOLERANCE_PERCENT = 0.30; // 30% relative error diff --git a/Test_Numerics/Distributions/Univariate/Test_Exponential.cs b/Test_Numerics/Distributions/Univariate/Test_Exponential.cs index c5ba0809..3113712f 100644 --- a/Test_Numerics/Distributions/Univariate/Test_Exponential.cs +++ b/Test_Numerics/Distributions/Univariate/Test_Exponential.cs @@ -1,4 +1,4 @@ -using System; +using System; using Microsoft.VisualStudio.TestTools.UnitTesting; using Numerics.Distributions; @@ -13,7 +13,7 @@ namespace Distributions.Univariate /// /// Haden Smith, USACE Risk Management Center, cole.h.smith@usace.army.mil /// Tiki Gonzalez, USACE Risk Management Center, julian.t.gonzalez@usace.army.mil - /// + /// /// /// /// References: @@ -22,7 +22,7 @@ namespace Distributions.Univariate /// /// /// - + [TestClass] public class Test_Exponential { @@ -35,7 +35,7 @@ public class Test_Exponential /// /// /// - /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. + /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. /// /// /// Example 6.1.1 page 132. @@ -81,7 +81,7 @@ public void Test_EXP_LMOM_Fit() /// /// /// - /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. + /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. /// /// /// Example 6.1.1 page 132. @@ -105,7 +105,7 @@ public void Test_EXP_MLE_Fit() /// /// /// - /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. + /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. /// /// /// Example 6.1.2 page 134. @@ -128,7 +128,7 @@ public void Test_EXP_Quantile() /// /// /// - /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. + /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. /// /// /// Example 6.1.3 page 138. diff --git a/Test_Numerics/Distributions/Univariate/Test_GammaDistribution.cs b/Test_Numerics/Distributions/Univariate/Test_GammaDistribution.cs index 13d6ce49..998cf0df 100644 --- a/Test_Numerics/Distributions/Univariate/Test_GammaDistribution.cs +++ b/Test_Numerics/Distributions/Univariate/Test_GammaDistribution.cs @@ -1,4 +1,4 @@ -using System; +using System; using Microsoft.VisualStudio.TestTools.UnitTesting; using Numerics.Distributions; @@ -13,7 +13,7 @@ namespace Distributions.Univariate /// /// Haden Smith, USACE Risk Management Center, cole.h.smith@usace.army.mil /// Tiki Gonzalez, USACE Risk Management Center, julian.t.gonzalez@usace.army.mil - /// + /// /// /// /// References: @@ -49,6 +49,9 @@ public void Test_GammaDist_MOM() Assert.IsLessThan(0.01d, (lambda - trueL) / trueL); } + /// + /// Verifies gamma distribution fitting with linear moments. + /// [TestMethod()] public void Test_GammaDist_LMOM_Fit() { @@ -74,7 +77,7 @@ public void Test_GammaDist_LMOM_Fit() /// /// /// - /// Reference: "The Gamma Family and Derived Distributions Applied in Hydrology", B. Bobee & F. Ashkar, Water Resources Publications, 1991. + /// Reference: "The Gamma Family and Derived Distributions Applied in Hydrology", B. Bobee & F. Ashkar, Water Resources Publications, 1991. /// /// /// Example 7.4 page 93. @@ -98,7 +101,7 @@ public void Test_GammaDist_MLE() /// /// /// - /// Reference: "The Gamma Family and Derived Distributions Applied in Hydrology", B. Bobee & F. Ashkar, Water Resources Publications, 1991. + /// Reference: "The Gamma Family and Derived Distributions Applied in Hydrology", B. Bobee & F. Ashkar, Water Resources Publications, 1991. /// /// /// Example 5.3 page 52. @@ -121,7 +124,7 @@ public void Test_GammaDist_Quantile() /// /// /// - /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. + /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. /// /// /// Example 6.1.3 page 138. @@ -319,7 +322,7 @@ public void Test_Maximum() /// /// /// - /// Reference: "The Gamma Family and Derived Distributions Applied in Hydrology", B. Bobee & F. Ashkar, Water Resources Publications, 1991. + /// Reference: "The Gamma Family and Derived Distributions Applied in Hydrology", B. Bobee & F. Ashkar, Water Resources Publications, 1991. /// /// /// Example 7.4 page 93. @@ -343,7 +346,7 @@ public void ValidateMLE_NR() /// /// /// - /// Reference: "The Gamma Family and Derived Distributions Applied in Hydrology", B. Bobee & F. Ashkar, Water Resources Publications, 1991. + /// Reference: "The Gamma Family and Derived Distributions Applied in Hydrology", B. Bobee & F. Ashkar, Water Resources Publications, 1991. /// /// /// Example 7.4 page 93. @@ -412,7 +415,7 @@ public void Test_InverseCDF() } /// - /// Validating Wilson Hilferty Inverse CDF function. + /// Validating Wilson Hilferty Inverse CDF function. /// /// [TestMethod()] @@ -426,5 +429,27 @@ public void ValidateWilsonHilfertyInverseCDF() Assert.AreEqual(0.3566877, G2.WilsonHilfertyInverseCDF(0.05),1e-04); Assert.AreEqual(0.5326, G2.WilsonHilfertyInverseCDF(0.10), 1e-04); } + + /// + /// Verifies the frequency-factor skew derivative at and around zero skew. + /// + [TestMethod] + public void Test_PartialKp_ZeroSkewLimitAndContinuity() + { + const double probability = 0.9d; + double z = Normal.StandardZ(probability); + double expected = (z * z - 1.0d) / 6.0d; + + Assert.AreEqual(expected, GammaDistribution.PartialKp(0.0d, probability), 1E-15); + Assert.AreEqual(expected, GammaDistribution.PartialKp(9.9E-5, probability), 1E-15); + Assert.AreEqual(expected, GammaDistribution.PartialKp(-9.9E-5, probability), 1E-15); + Assert.AreEqual(expected, GammaDistribution.PartialKp(1.01E-4, probability), 1E-4); + Assert.AreEqual(expected, GammaDistribution.PartialKp(-1.01E-4, probability), 1E-4); + + const double h = 1E-3; + double finiteDifference = (GammaDistribution.FrequencyFactorKp(h, probability) - + GammaDistribution.FrequencyFactorKp(-h, probability)) / (2.0d * h); + Assert.AreEqual(expected, finiteDifference, 1E-4); + } } } diff --git a/Test_Numerics/Distributions/Univariate/Test_GeneralizedBeta.cs b/Test_Numerics/Distributions/Univariate/Test_GeneralizedBeta.cs index 1ff56d36..bec72ce8 100644 --- a/Test_Numerics/Distributions/Univariate/Test_GeneralizedBeta.cs +++ b/Test_Numerics/Distributions/Univariate/Test_GeneralizedBeta.cs @@ -38,7 +38,7 @@ public void Test_GenBeta() double true_mean = 0.21105527638190955d; double true_median = 0.11577706212908731d; - double true_mode = 57.999999999999957d; + double true_mode = 0.0d; double true_var = 0.055689279830523512d; double true_pdf = 0.94644031936694828d; double true_cdf = 0.69358638272337991d; @@ -332,5 +332,27 @@ public void Test_Mode_UniformCase() var b4 = new GeneralizedBeta(5, 2, 0, 1); Assert.AreEqual(0.8, b4.Mode, 1E-10); } + + /// + /// Verifies boundary and legacy non-unique mode conventions. + /// + [TestMethod] + public void Test_Mode_BoundaryShapes() + { + var decreasing = new GeneralizedBeta(0.42d, 1.57d, -10.0d, 20.0d); + Assert.AreEqual(-10.0d, decreasing.Mode); + + var increasing = new GeneralizedBeta(1.57d, 0.42d, -10.0d, 20.0d); + Assert.AreEqual(20.0d, increasing.Mode); + + var lowerBoundary = new GeneralizedBeta(0.5d, 1.0d, -10.0d, 20.0d); + Assert.AreEqual(-10.0d, lowerBoundary.Mode); + + var upperBoundary = new GeneralizedBeta(1.0d, 0.5d, -10.0d, 20.0d); + Assert.AreEqual(20.0d, upperBoundary.Mode); + + var uShaped = new GeneralizedBeta(0.5d, 0.5d, -10.0d, 20.0d); + Assert.AreEqual(5.0d, uShaped.Mode); + } } } diff --git a/Test_Numerics/Distributions/Univariate/Test_GeneralizedExtremeValue.cs b/Test_Numerics/Distributions/Univariate/Test_GeneralizedExtremeValue.cs index aa0e5a70..c2f161da 100644 --- a/Test_Numerics/Distributions/Univariate/Test_GeneralizedExtremeValue.cs +++ b/Test_Numerics/Distributions/Univariate/Test_GeneralizedExtremeValue.cs @@ -1,4 +1,4 @@ -using System; +using System; using System.Linq; using Microsoft.VisualStudio.TestTools.UnitTesting; using Numerics; @@ -15,7 +15,7 @@ namespace Distributions.Univariate /// /// Haden Smith, USACE Risk Management Center, cole.h.smith@usace.army.mil /// Tiki Gonzalez, USACE Risk Management Center, julian.t.gonzalez@usace.army.mil - /// + /// /// /// /// References: @@ -37,7 +37,7 @@ public class Test_GeneralizedExtremeValue /// /// /// - /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. + /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. /// /// /// Example 7.1.1 page 218. @@ -64,7 +64,7 @@ public void Test_GEV_MOM_Fit() /// /// /// - /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. + /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. /// /// /// Example 7.1.1 page 218. @@ -97,7 +97,7 @@ public void Test_GEV_LMOM_Fit() /// /// /// - /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. + /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. /// /// /// Example 7.1.1 page 219. @@ -124,7 +124,7 @@ public void Test_GEV_MLE_Fit() /// /// /// - /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. + /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. /// /// /// Example 7.1.2 page 221. @@ -147,7 +147,7 @@ public void Test_GEV_Quantile() /// /// /// - /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. + /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. /// /// /// Example 7.1.3 page 226. diff --git a/Test_Numerics/Distributions/Univariate/Test_GeneralizedLogistic.cs b/Test_Numerics/Distributions/Univariate/Test_GeneralizedLogistic.cs index c2beca91..ea7f213d 100644 --- a/Test_Numerics/Distributions/Univariate/Test_GeneralizedLogistic.cs +++ b/Test_Numerics/Distributions/Univariate/Test_GeneralizedLogistic.cs @@ -1,4 +1,4 @@ -using Microsoft.VisualStudio.TestTools.UnitTesting; +using Microsoft.VisualStudio.TestTools.UnitTesting; using Numerics.Distributions; namespace Distributions.Univariate @@ -12,7 +12,7 @@ namespace Distributions.Univariate /// /// Haden Smith, USACE Risk Management Center, cole.h.smith@usace.army.mil /// Tiki Gonzalez, USACE Risk Management Center, julian.t.gonzalez@usace.army.mil - /// + /// /// /// /// References: @@ -34,7 +34,7 @@ public class Test_GeneralizedLogistic /// /// /// - /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. + /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. /// /// /// Example 9.2.1 page 311. @@ -64,7 +64,7 @@ public void Test_GLO_MOM_Fit() /// /// /// - /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. + /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. /// /// /// Example 9.2.1 page 311. @@ -97,7 +97,7 @@ public void Test_GLO_LMOM_Fit() /// /// /// - /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. + /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. /// /// /// Example 9.1.1 page 313. @@ -124,7 +124,7 @@ public void Test_GLO_MLE_Fit() /// /// /// - /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. + /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. /// /// /// Example 9.2.2 page 315. @@ -263,7 +263,7 @@ public void Test_Mode() var l2 = new GeneralizedLogistic(100, 10, 1); Assert.AreEqual(95, l2.Mode); } - + /// /// Checking Standard deviation. /// @@ -329,7 +329,7 @@ public void Test_Minimum() /// Testing maximum function. /// [TestMethod()] - public void Test_Maximum() + public void Test_Maximum() { var l = new GeneralizedLogistic(); Assert.AreEqual(double.PositiveInfinity, l.Maximum); @@ -382,5 +382,42 @@ public void Test_InverseCDF() Assert.AreEqual(100, l2.InverseCDF(0.5)); Assert.AreEqual(105.714285, l2.InverseCDF(0.7), 1e-04); } + + /// + /// Verifies the Logistic L-moment limit and stable near-zero round trips. + /// + [TestMethod] + public void Test_LinearMoments_ZeroAndNearZeroKappa() + { + var distribution = new GeneralizedLogistic(); + double[] logistic = distribution.LinearMomentsFromParameters([2.0d, 3.0d, 0.0d]); + + Assert.AreEqual(2.0d, logistic[0]); + Assert.AreEqual(3.0d, logistic[1]); + Assert.AreEqual(0.0d, logistic[2]); + Assert.AreEqual(1.0d / 6.0d, logistic[3], 1E-15); + + double[] recoveredLogistic = distribution.ParametersFromLinearMoments(logistic); + Assert.AreEqual(2.0d, recoveredLogistic[0]); + Assert.AreEqual(3.0d, recoveredLogistic[1]); + Assert.AreEqual(0.0d, recoveredLogistic[2]); + + foreach (double kappa in new[] { -1E-8, 1E-8, -5E-5, 5E-5 }) + { + double[] parameters = [2.0d, 3.0d, kappa]; + double[] moments = distribution.LinearMomentsFromParameters(parameters); + double[] recovered = distribution.ParametersFromLinearMoments(moments); + + for (int i = 0; i < moments.Length; i++) + { + Assert.IsFalse(double.IsNaN(moments[i])); + Assert.IsFalse(double.IsInfinity(moments[i])); + } + + Assert.AreEqual(parameters[0], recovered[0], 1E-12); + Assert.AreEqual(parameters[1], recovered[1], 1E-12); + Assert.AreEqual(parameters[2], recovered[2], 1E-15); + } + } } } diff --git a/Test_Numerics/Distributions/Univariate/Test_GeneralizedPareto.cs b/Test_Numerics/Distributions/Univariate/Test_GeneralizedPareto.cs index 5687846c..c1621096 100644 --- a/Test_Numerics/Distributions/Univariate/Test_GeneralizedPareto.cs +++ b/Test_Numerics/Distributions/Univariate/Test_GeneralizedPareto.cs @@ -1,4 +1,4 @@ -using System; +using System; using System.Diagnostics; using Microsoft.VisualStudio.TestTools.UnitTesting; using Numerics.Data.Statistics; @@ -15,7 +15,7 @@ namespace Distributions.Univariate /// /// Haden Smith, USACE Risk Management Center, cole.h.smith@usace.army.mil /// Tiki Gonzalez, USACE Risk Management Center, julian.t.gonzalez@usace.army.mil - /// + /// /// /// /// References: @@ -37,7 +37,7 @@ public class Test_GeneralizedPareto /// /// /// - /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. + /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. /// /// /// Example 8.3.1 page 279. @@ -89,7 +89,7 @@ public void Test_GPA_LMOM_Fit() /// /// /// - /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. + /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. /// /// /// Example 8.3.1 page 279. @@ -116,7 +116,7 @@ public void Test_GPA_ModMOM_Fit() /// /// /// - /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. + /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. /// /// /// Example 8.3.1 page 279. @@ -143,7 +143,7 @@ public void Test_GPA_MLE_Fit() /// /// /// - /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. + /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. /// /// /// Example 8.3.2 page 283. @@ -184,7 +184,7 @@ public void Test_GPA_Partials() /// /// /// - /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. + /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. /// /// /// Example 8.3.3 page 286. @@ -326,7 +326,7 @@ public void Test_StandardDeviation() var GPA3 = new GeneralizedPareto(100, 10, 1); Assert.AreEqual(double.NaN, GPA3.StandardDeviation); } - + /// /// Testing the skew function. /// @@ -373,7 +373,7 @@ public void Test_Minimum() /// Testing maximum function. /// [TestMethod()] - public void Test_Maximum() + public void Test_Maximum() { var GPA = new GeneralizedPareto(); Assert.AreEqual(double.PositiveInfinity, GPA.Maximum); diff --git a/Test_Numerics/Distributions/Univariate/Test_Gumbel.cs b/Test_Numerics/Distributions/Univariate/Test_Gumbel.cs index 55766142..0c768baf 100644 --- a/Test_Numerics/Distributions/Univariate/Test_Gumbel.cs +++ b/Test_Numerics/Distributions/Univariate/Test_Gumbel.cs @@ -1,4 +1,4 @@ -using System; +using System; using Microsoft.VisualStudio.TestTools.UnitTesting; using Numerics.Distributions; @@ -13,7 +13,7 @@ namespace Distributions.Univariate /// /// Haden Smith, USACE Risk Management Center, cole.h.smith@usace.army.mil /// Tiki Gonzalez, USACE Risk Management Center, julian.t.gonzalez@usace.army.mil - /// + /// /// /// /// References: @@ -35,7 +35,7 @@ public class Test_Gumbel /// /// /// - /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. + /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. /// /// /// Example 7.2.1 page 234. @@ -81,7 +81,7 @@ public void Test_GUM_LMOM_Fit() /// /// /// - /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. + /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. /// /// /// Example 7.2.1 page 234. @@ -105,7 +105,7 @@ public void Test_GUM_MLE_Fit() /// /// /// - /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. + /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. /// /// /// Example 7.7.2 page 237. @@ -129,7 +129,7 @@ public void Test_Gumbel_Quantile() /// /// /// - /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. + /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. /// /// /// Example 7.2.3 page 240. diff --git a/Test_Numerics/Distributions/Univariate/Test_KappaFour.cs b/Test_Numerics/Distributions/Univariate/Test_KappaFour.cs index 2ef0162c..f48594a8 100644 --- a/Test_Numerics/Distributions/Univariate/Test_KappaFour.cs +++ b/Test_Numerics/Distributions/Univariate/Test_KappaFour.cs @@ -1,4 +1,4 @@ -using Microsoft.VisualStudio.TestTools.UnitTesting; +using Microsoft.VisualStudio.TestTools.UnitTesting; using Numerics.Distributions; using Numerics.Mathematics; @@ -13,7 +13,7 @@ namespace Distributions.Univariate /// /// Haden Smith, USACE Risk Management Center, cole.h.smith@usace.army.mil /// Tiki Gonzalez, USACE Risk Management Center, julian.t.gonzalez@usace.army.mil - /// + /// /// /// /// References: @@ -26,6 +26,9 @@ namespace Distributions.Univariate public class Test_KappaFour { + /// + /// Verifies Kappa Four fitting with linear moments. + /// [TestMethod] public void Test_K4_LMOM() { @@ -55,7 +58,7 @@ public void Test_K4_LMOM() } /// - /// Verification using R lmom package. + /// Verification using R lmom package. /// [TestMethod] public void Test_K4_CDF() @@ -77,6 +80,9 @@ public void Test_K4_CDF() } + /// + /// Verifies Kappa Four distribution values. + /// [TestMethod] public void Test_K4_Dist() { @@ -98,10 +104,13 @@ public void Test_K4_Dist() } + /// + /// Verifies Kappa Four partial derivative calculations. + /// [TestMethod] public void Test_K4_PartialDerivatives() { - + // Air quality - wind data from R var data = new double[] { 7.4, 8, 12.6, 11.5, 14.3, 14.9, 8.6, 13.8, 20.1, 8.6, 6.9, 9.7, 9.2, 10.9, 13.2, 11.5, 12, 18.4, 11.5, 9.7, 9.7, 16.6, 9.7, 12, 16.6, 14.9, 8, 12, 14.9, 5.7, 7.4, 8.6, 9.7, 16.1, 9.2, 8.6, 14.3, 9.7, 6.9, 13.8, 11.5, 10.9, 9.2, 8, 13.8, 11.5, 14.9, 20.7, 9.2, 11.5, 10.3, 6.3, 1.7, 4.6, 6.3, 8, 8, 10.3, 11.5, 14.9, 8, 4.1, 9.2, 9.2, 10.9, 4.6, 10.9, 5.1, 6.3, 5.7, 7.4, 8.6, 14.3, 14.9, 14.9, 14.3, 6.9, 10.3, 6.3, 5.1, 11.5, 6.9, 9.7, 11.5, 8.6, 8, 8.6, 12, 7.4, 7.4, 7.4, 9.2, 6.9, 13.8, 7.4, 6.9, 7.4, 4.6, 4, 10.3, 8, 8.6, 11.5, 11.5, 11.5, 9.7, 11.5, 10.3, 6.3, 7.4, 10.9, 10.3, 15.5, 14.3, 12.6, 9.7, 3.4, 8, 5.7, 9.7, 2.3, 6.3, 6.3, 6.9, 5.1, 2.8, 4.6, 7.4, 15.5, 10.9, 10.3, 10.9, 9.7, 14.9, 15.5, 6.3, 10.9, 11.5, 6.9, 13.8, 10.3, 10.3, 8, 12.6, 9.2, 10.3, 10.3, 16.6, 6.9, 13.2, 14.3, 8, 11.5 }; var kappa4 = new KappaFour(); @@ -109,7 +118,7 @@ public void Test_K4_PartialDerivatives() double p = 0.999; var pd1 = kappa4.QuantileGradient(p); - var pd2 = NumericalDerivative.Gradient(x => + var pd2 = NumericalDerivative.Gradient(x => { var K4 = new KappaFour(); K4.SetParameters(x); diff --git a/Test_Numerics/Distributions/Univariate/Test_LnNormal.cs b/Test_Numerics/Distributions/Univariate/Test_LnNormal.cs index c8c60bef..f1b7d0bd 100644 --- a/Test_Numerics/Distributions/Univariate/Test_LnNormal.cs +++ b/Test_Numerics/Distributions/Univariate/Test_LnNormal.cs @@ -1,4 +1,4 @@ -using System; +using System; using Microsoft.VisualStudio.TestTools.UnitTesting; using Numerics.Distributions; @@ -13,7 +13,7 @@ namespace Distributions.Univariate /// /// Haden Smith, USACE Risk Management Center, cole.h.smith@usace.army.mil /// Tiki Gonzalez, USACE Risk Management Center, julian.t.gonzalez@usace.army.mil - /// + /// /// /// /// References: @@ -35,7 +35,7 @@ public class Test_LnNormal /// /// /// - /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. + /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. /// /// /// Example 5.2.1 page 99. @@ -60,7 +60,7 @@ public void Test_LnNormal_MOM_Fit() /// /// /// - /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. + /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. /// /// /// Example 5.2.1 page 99. @@ -84,7 +84,7 @@ public void Test_LnNormal_MLE_Fit() /// /// /// - /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. + /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. /// /// /// Example 5.2.2 page 102. @@ -107,7 +107,7 @@ public void Test_LnNormal_Quantile() /// /// /// - /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. + /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. /// /// /// Example 5.1.3 page 94. diff --git a/Test_Numerics/Distributions/Univariate/Test_LogNormal.cs b/Test_Numerics/Distributions/Univariate/Test_LogNormal.cs index 3857363d..a33b0636 100644 --- a/Test_Numerics/Distributions/Univariate/Test_LogNormal.cs +++ b/Test_Numerics/Distributions/Univariate/Test_LogNormal.cs @@ -1,4 +1,4 @@ -using System; +using System; using Microsoft.VisualStudio.TestTools.UnitTesting; using Numerics.Distributions; @@ -13,7 +13,7 @@ namespace Distributions.Univariate /// /// Haden Smith, USACE Risk Management Center, cole.h.smith@usace.army.mil /// Tiki Gonzalez, USACE Risk Management Center, julian.t.gonzalez@usace.army.mil - /// + /// /// /// /// References: @@ -35,7 +35,7 @@ public class Test_LogNormal /// /// /// - /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. + /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. /// /// /// Example 5.2.1 page 99. @@ -82,7 +82,7 @@ public void Test_LogNormal_LMOM_Fit() /// /// /// - /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. + /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. /// /// /// Example 5.2.1 page 99. @@ -108,7 +108,7 @@ public void Test_LogNormal_MLE_Fit() /// /// /// - /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. + /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. /// /// /// Example 5.2.2 page 102. @@ -131,7 +131,7 @@ public void Test_LogNormal_Quantile() /// /// /// - /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. + /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. /// /// /// Example 5.1.3 page 94. @@ -171,7 +171,7 @@ public void Test_InvalidParameters() { var LogN = new LogNormal(double.NaN,double.NaN); Assert.IsFalse(LogN.ParametersValid); - + var LogN2 = new LogNormal(double.PositiveInfinity,double.PositiveInfinity); Assert.IsFalse(LogN2.ParametersValid); diff --git a/Test_Numerics/Distributions/Univariate/Test_LogPearsonTypeIII.cs b/Test_Numerics/Distributions/Univariate/Test_LogPearsonTypeIII.cs index fef5c754..a59e223b 100644 --- a/Test_Numerics/Distributions/Univariate/Test_LogPearsonTypeIII.cs +++ b/Test_Numerics/Distributions/Univariate/Test_LogPearsonTypeIII.cs @@ -1,4 +1,4 @@ -using System; +using System; using Microsoft.VisualStudio.TestTools.UnitTesting; using Numerics.Distributions; @@ -13,7 +13,7 @@ namespace Distributions.Univariate /// /// Haden Smith, USACE Risk Management Center, cole.h.smith@usace.army.mil /// Tiki Gonzalez, USACE Risk Management Center, julian.t.gonzalez@usace.army.mil - /// + /// /// /// /// References: @@ -35,7 +35,7 @@ public class Test_LogPearsonTypeIII /// /// /// - /// Reference: "The Gamma Family and Derived Distributions Applied in Hydrology", B. Bobee & F. Ashkar, Water Resources Publications, 1991. + /// Reference: "The Gamma Family and Derived Distributions Applied in Hydrology", B. Bobee & F. Ashkar, Water Resources Publications, 1991. /// /// /// Example 7.4 page 93. @@ -80,7 +80,7 @@ public void Test_LP3_IndirectMOM() /// /// /// - /// Reference: "The Gamma Family and Derived Distributions Applied in Hydrology", B. Bobee & F. Ashkar, Water Resources Publications, 1991. + /// Reference: "The Gamma Family and Derived Distributions Applied in Hydrology", B. Bobee & F. Ashkar, Water Resources Publications, 1991. /// /// /// Example 7.4 page 93. @@ -125,7 +125,7 @@ public void Test_LP3_MLE() /// /// /// - /// Reference: "The Gamma Family and Derived Distributions Applied in Hydrology", B. Bobee & F. Ashkar, Water Resources Publications, 1991. + /// Reference: "The Gamma Family and Derived Distributions Applied in Hydrology", B. Bobee & F. Ashkar, Water Resources Publications, 1991. /// /// /// Example 7.1 page 87. @@ -148,10 +148,10 @@ public void Test_LP3_Quantile() /// /// /// - /// Reference: "The Gamma Family and Derived Distributions Applied in Hydrology", B. Bobee & F. Ashkar, Water Resources Publications, 1991. + /// Reference: "The Gamma Family and Derived Distributions Applied in Hydrology", B. Bobee & F. Ashkar, Water Resources Publications, 1991. /// /// - /// Example 7.1 & 6.4 page 87-93. + /// Example 7.1 & 6.4 page 87-93. /// /// [TestMethod()] @@ -269,7 +269,7 @@ public void Test_Minimum() /// Testing maximum. /// [TestMethod()] - public void Test_Maximum() + public void Test_Maximum() { var LP3 = new LogPearsonTypeIII(); Assert.AreEqual(double.PositiveInfinity,LP3.Maximum ); @@ -284,7 +284,7 @@ public void Test_Maximum() /// /// Testing PDF method. /// - [TestMethod()] + [TestMethod()] public void Test_PDF() { var LP3 = new LogPearsonTypeIII(); @@ -326,5 +326,43 @@ public void ValidateWilsonHilfertyInverseCDF() Assert.AreEqual(double.PositiveInfinity, LP3.WilsonHilfertyInverseCDF(1)); Assert.AreEqual(747.01005, LP3.WilsonHilfertyInverseCDF(0.4), 1e-05); } + + /// + /// Verifies signed, finite Log-Pearson III L-moments and their Normal limit. + /// + [TestMethod] + public void Test_LinearMoments_SignedSmallAndZeroSkew() + { + var distribution = new LogPearsonTypeIII(); + double[] normalMoments = distribution.LinearMomentsFromParameters([2.0d, 0.3d, 0.0d]); + Assert.AreEqual(2.0d, normalMoments[0]); + Assert.AreEqual(0.3d / Math.Sqrt(Math.PI), normalMoments[1], 1E-14); + Assert.AreEqual(0.0d, normalMoments[2]); + Assert.AreEqual(0.12260172d, normalMoments[3], 1E-14); + + double[] normalParameters = distribution.ParametersFromLinearMoments(normalMoments); + Assert.AreEqual(2.0d, normalParameters[0]); + Assert.AreEqual(0.3d, normalParameters[1], 1E-14); + Assert.AreEqual(0.0d, normalParameters[2]); + + double[] positive = distribution.LinearMomentsFromParameters([2.0d, 0.3d, 0.1d]); + double[] negative = distribution.LinearMomentsFromParameters([2.0d, 0.3d, -0.1d]); + Assert.AreEqual(2.0d, positive[0]); + Assert.AreEqual(positive[0], negative[0]); + Assert.AreEqual(positive[1], negative[1], 1E-14); + Assert.AreEqual(positive[2], -negative[2], 1E-14); + Assert.IsGreaterThan(0.0d, positive[2]); + + foreach (double moment in positive) + { + Assert.IsFalse(double.IsNaN(moment)); + Assert.IsFalse(double.IsInfinity(moment)); + } + + double[] recovered = distribution.ParametersFromLinearMoments(negative); + Assert.AreEqual(2.0d, recovered[0]); + Assert.AreEqual(0.3d, recovered[1], 1E-5); + Assert.AreEqual(-0.1d, recovered[2], 1E-4); + } } } diff --git a/Test_Numerics/Distributions/Univariate/Test_Logistic.cs b/Test_Numerics/Distributions/Univariate/Test_Logistic.cs index 456a503f..1d8a02ef 100644 --- a/Test_Numerics/Distributions/Univariate/Test_Logistic.cs +++ b/Test_Numerics/Distributions/Univariate/Test_Logistic.cs @@ -1,4 +1,4 @@ -using System; +using System; using Microsoft.VisualStudio.TestTools.UnitTesting; using Numerics.Distributions; @@ -13,7 +13,7 @@ namespace Distributions.Univariate /// /// Haden Smith, USACE Risk Management Center, cole.h.smith@usace.army.mil /// Tiki Gonzalez, USACE Risk Management Center, julian.t.gonzalez@usace.army.mil - /// + /// /// /// /// References: @@ -35,7 +35,7 @@ public class Test_Logistic /// /// /// - /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. + /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. /// /// /// Example 9.1.1 page 295. @@ -59,7 +59,7 @@ public void Test_Logistic_MOM_Fit() /// /// /// - /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. + /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. /// /// /// Example 9.1.1 page 295. @@ -83,7 +83,7 @@ public void Test_Logistic_MLE_Fit() /// /// /// - /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. + /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. /// /// /// Example 9.1.2 page 297. @@ -106,7 +106,7 @@ public void Test_Logistic_Quantile() /// /// /// - /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. + /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. /// /// /// Example 9.1.3 page 300. diff --git a/Test_Numerics/Distributions/Univariate/Test_NoncentralT.cs b/Test_Numerics/Distributions/Univariate/Test_NoncentralT.cs index 61507a69..e6f5a94c 100644 --- a/Test_Numerics/Distributions/Univariate/Test_NoncentralT.cs +++ b/Test_Numerics/Distributions/Univariate/Test_NoncentralT.cs @@ -77,6 +77,49 @@ public void Test_NoncentralT_InverseCDF() } + /// + /// Verifies that zero noncentrality reduces to the central Student-t distribution across both tails. + /// + [TestMethod] + public void Test_NoncentralT_ZeroNoncentralityMatchesStudentT() + { + var noncentral = new NoncentralT(12d, 0d); + var central = new StudentT(0d, 1d, 12d); + foreach (double x in new[] { -3d, -1d, 0d, 1d, 3d }) + { + Assert.AreEqual(central.CDF(x), noncentral.CDF(x), 1E-8, $"CDF mismatch at x={x}."); + Assert.AreEqual(central.PDF(x), noncentral.PDF(x), 1E-7, $"PDF mismatch at x={x}."); + } + } + + /// + /// Verifies negative-tail symmetry, the exact value at zero, and inverse-CDF round trips. + /// + [TestMethod] + public void Test_NoncentralT_NegativeTailSymmetryAndRoundTrips() + { + const double degreesOfFreedom = 9d; + const double noncentrality = 1.75d; + var positive = new NoncentralT(degreesOfFreedom, noncentrality); + var reflected = new NoncentralT(degreesOfFreedom, -noncentrality); + + foreach (double x in new[] { -4d, -2d, -0.5d, 0d, 1.5d }) + { + Assert.AreEqual( + positive.CDF(x), + 1d - reflected.CDF(-x), + 1E-8, + $"Reflection identity failed at x={x}."); + } + Assert.AreEqual(Normal.StandardCDF(-noncentrality), positive.CDF(0d), 1E-12); + + foreach (double probability in new[] { 0.05d, 0.25d, 0.75d, 0.95d }) + { + double quantile = positive.InverseCDF(probability); + Assert.AreEqual(probability, positive.CDF(quantile), 1E-6, $"Round trip failed at p={probability}."); + } + } + /// /// Verifying input parameters can create Noncentral T distribution. /// diff --git a/Test_Numerics/Distributions/Univariate/Test_Normal.cs b/Test_Numerics/Distributions/Univariate/Test_Normal.cs index d28289f4..8cf140f6 100644 --- a/Test_Numerics/Distributions/Univariate/Test_Normal.cs +++ b/Test_Numerics/Distributions/Univariate/Test_Normal.cs @@ -1,4 +1,4 @@ -using System; +using System; using Microsoft.VisualStudio.TestTools.UnitTesting; using Numerics.Distributions; @@ -13,7 +13,7 @@ namespace Distributions.Univariate /// /// Haden Smith, USACE Risk Management Center, cole.h.smith@usace.army.mil /// Tiki Gonzalez, USACE Risk Management Center, julian.t.gonzalez@usace.army.mil - /// + /// /// /// /// References: @@ -35,7 +35,7 @@ public class Test_Normal /// /// /// - /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. + /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. /// /// /// Example 5.1.1 page 88. @@ -106,7 +106,7 @@ public void Test_Normal_MLE_Fit() /// /// /// - /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. + /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. /// /// /// Example 5.1.2 page 91. @@ -129,7 +129,7 @@ public void Test_Normal_Quantile() /// /// /// - /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. + /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. /// /// /// Example 5.1.3 page 94. diff --git a/Test_Numerics/Distributions/Univariate/Test_PearsonTypeIII.cs b/Test_Numerics/Distributions/Univariate/Test_PearsonTypeIII.cs index ad2da5c6..c11f3d58 100644 --- a/Test_Numerics/Distributions/Univariate/Test_PearsonTypeIII.cs +++ b/Test_Numerics/Distributions/Univariate/Test_PearsonTypeIII.cs @@ -1,4 +1,4 @@ -using System; +using System; using Microsoft.VisualStudio.TestTools.UnitTesting; using Numerics.Distributions; @@ -13,7 +13,7 @@ namespace Distributions.Univariate /// /// Haden Smith, USACE Risk Management Center, cole.h.smith@usace.army.mil /// Tiki Gonzalez, USACE Risk Management Center, julian.t.gonzalez@usace.army.mil - /// + /// /// /// /// References: @@ -35,7 +35,7 @@ public class Test_PearsonTypeIII /// /// /// - /// Reference: "The Gamma Family and Derived Distributions Applied in Hydrology", B. Bobee & F. Ashkar, Water Resources Publications, 1991. + /// Reference: "The Gamma Family and Derived Distributions Applied in Hydrology", B. Bobee & F. Ashkar, Water Resources Publications, 1991. /// /// /// Example 6.3 page 70. @@ -105,7 +105,7 @@ public void Test_P3_LMOM_Fit() /// /// /// - /// Reference: "The Gamma Family and Derived Distributions Applied in Hydrology", B. Bobee & F. Ashkar, Water Resources Publications, 1991. + /// Reference: "The Gamma Family and Derived Distributions Applied in Hydrology", B. Bobee & F. Ashkar, Water Resources Publications, 1991. /// /// /// Example 6.1 page 64. @@ -150,7 +150,7 @@ public void Test_P3_MLE() /// /// /// - /// Reference: "The Gamma Family and Derived Distributions Applied in Hydrology", B. Bobee & F. Ashkar, Water Resources Publications, 1991. + /// Reference: "The Gamma Family and Derived Distributions Applied in Hydrology", B. Bobee & F. Ashkar, Water Resources Publications, 1991. /// /// /// Example 6.1 page 64. @@ -170,10 +170,10 @@ public void Test_P3_Quantile() /// /// /// - /// Reference: "The Gamma Family and Derived Distributions Applied in Hydrology", B. Bobee & F. Ashkar, Water Resources Publications, 1991. + /// Reference: "The Gamma Family and Derived Distributions Applied in Hydrology", B. Bobee & F. Ashkar, Water Resources Publications, 1991. /// /// - /// Example 6.1 & 6.3 page 64-70. + /// Example 6.1 & 6.3 page 64-70. /// /// [TestMethod()] @@ -347,7 +347,7 @@ public void Test_Minimum() /// Testing maximum function. /// [TestMethod()] - public void Test_Maximum() + public void Test_Maximum() { var P3 = new PearsonTypeIII(); Assert.AreEqual(double.PositiveInfinity, P3.Maximum); @@ -356,6 +356,48 @@ public void Test_Maximum() Assert.AreEqual(double.PositiveInfinity, P3ii.Maximum); } + /// + /// Verifies signed, finite Pearson III L-moments and their Normal limit. + /// + [TestMethod] + public void Test_LinearMoments_SignedSmallAndZeroSkew() + { + var distribution = new PearsonTypeIII(); + double[] normalMoments = distribution.LinearMomentsFromParameters([10.0d, 2.0d, 0.0d]); + Assert.AreEqual(10.0d, normalMoments[0]); + Assert.AreEqual(2.0d / Math.Sqrt(Math.PI), normalMoments[1], 1E-14); + Assert.AreEqual(0.0d, normalMoments[2]); + Assert.AreEqual(0.12260172d, normalMoments[3], 1E-14); + + double[] normalParameters = distribution.ParametersFromLinearMoments(normalMoments); + Assert.AreEqual(10.0d, normalParameters[0]); + Assert.AreEqual(2.0d, normalParameters[1], 1E-14); + Assert.AreEqual(0.0d, normalParameters[2]); + + double[] positive = distribution.LinearMomentsFromParameters([10.0d, 2.0d, 0.1d]); + double[] negative = distribution.LinearMomentsFromParameters([10.0d, 2.0d, -0.1d]); + Assert.AreEqual(10.0d, positive[0]); + Assert.AreEqual(positive[0], negative[0]); + Assert.AreEqual(positive[1], negative[1], 1E-14); + Assert.AreEqual(positive[2], -negative[2], 1E-14); + Assert.IsGreaterThan(0.0d, positive[2]); + + foreach (double moment in positive) + { + Assert.IsFalse(double.IsNaN(moment)); + Assert.IsFalse(double.IsInfinity(moment)); + } + + double[] recovered = distribution.ParametersFromLinearMoments(negative); + Assert.AreEqual(10.0d, recovered[0]); + Assert.AreEqual(2.0d, recovered[1], 1E-4); + Assert.AreEqual(-0.1d, recovered[2], 1E-4); + + double expectedDerivative = 2.0d * GammaDistribution.PartialKp(0.0d, 0.9d); + var zeroSkew = new PearsonTypeIII(10.0d, 2.0d, 0.0d); + Assert.AreEqual(expectedDerivative, zeroSkew.QuantileGradientForMoments(0.9d)[2], 1E-14); + } + } diff --git a/Test_Numerics/Distributions/Univariate/Test_Pert.cs b/Test_Numerics/Distributions/Univariate/Test_Pert.cs index d97feb20..4394c278 100644 --- a/Test_Numerics/Distributions/Univariate/Test_Pert.cs +++ b/Test_Numerics/Distributions/Univariate/Test_Pert.cs @@ -255,5 +255,27 @@ public void Test_MinMax() Assert.AreEqual(1, p2.Minimum); Assert.AreEqual(2, p2.Maximum); } + + /// + /// Preserves the TotalRisk convention that equal PERT inputs represent a fixed value. + /// + [TestMethod] + public void Test_Pert_EqualInputsRemainFixedValue() + { + const double value = 42.0d; + var distribution = new Pert(value, value, value); + + Assert.IsTrue(distribution.ParametersValid); + Assert.AreEqual(value, distribution.Mean); + Assert.AreEqual(value, distribution.Median); + Assert.AreEqual(value, distribution.Mode); + Assert.AreEqual(0.0d, distribution.StandardDeviation); + + double[] probabilities = [0.0d, 0.01d, 0.5d, 0.99d, 1.0d]; + foreach (double probability in probabilities) + { + Assert.AreEqual(value, distribution.InverseCDF(probability)); + } + } } } diff --git a/Test_Numerics/Distributions/Univariate/Test_StudentT.cs b/Test_Numerics/Distributions/Univariate/Test_StudentT.cs index eb35d32e..1df25c08 100644 --- a/Test_Numerics/Distributions/Univariate/Test_StudentT.cs +++ b/Test_Numerics/Distributions/Univariate/Test_StudentT.cs @@ -36,7 +36,7 @@ public void Test_StudentT_PDF() Assert.AreEqual(result, pdf, 1E-10); t = new StudentT(2.5d, 0.5d, 4d); pdf = t.PDF(1.4d); - result = 0.0516476521260042d; + result = 0.1032953042520084d; Assert.AreEqual(result, pdf, 1E-10); } @@ -230,5 +230,48 @@ public void Test_MinMax() Assert.AreEqual(double.NegativeInfinity, t.Minimum); Assert.AreEqual(double.PositiveInfinity, t.Maximum); } + + /// + /// Verifies that location-scale Student-t densities include the scale Jacobian. + /// + [TestMethod] + public void Test_StudentT_PDF_LocationScaleIdentityAndNormalization() + { + var standard = new StudentT(0.0d, 1.0d, 4.0d); + var scaled = new StudentT(2.5d, 0.5d, 4.0d); + const double x = 1.4d; + double z = (x - scaled.Mu) / scaled.Sigma; + + Assert.AreEqual(standard.PDF(z) / scaled.Sigma, scaled.PDF(x), 1E-14); + + const int intervals = 100000; + double lower = scaled.Mu - 100.0d * scaled.Sigma; + double upper = scaled.Mu + 100.0d * scaled.Sigma; + double width = (upper - lower) / intervals; + double integral = 0.5d * (scaled.PDF(lower) + scaled.PDF(upper)); + for (int i = 1; i < intervals; i++) + { + integral += scaled.PDF(lower + i * width); + } + + Assert.AreEqual(1.0d, integral * width, 1E-6); + } + + /// + /// Verifies that extreme finite Student-t quantiles retain location and scale. + /// + [TestMethod] + public void Test_StudentT_InverseCDF_ExtremeTailRetainsLocationAndScale() + { + const double probability = 1E-155; + var distribution = new StudentT(2.0d, 3.0d, 1.0d); + + double actual = distribution.InverseCDF(probability); + double expected = 2.0d - 3.0d / (Math.PI * probability); + + Assert.IsFalse(double.IsInfinity(actual)); + Assert.AreEqual(expected, actual, Math.Abs(expected) * 1E-12); + Assert.IsLessThan(distribution.InverseCDF(1E-150), actual); + } } } diff --git a/Test_Numerics/Distributions/Univariate/Test_UnivariateDistributionFactory.cs b/Test_Numerics/Distributions/Univariate/Test_UnivariateDistributionFactory.cs new file mode 100644 index 00000000..7bffc26c --- /dev/null +++ b/Test_Numerics/Distributions/Univariate/Test_UnivariateDistributionFactory.cs @@ -0,0 +1,71 @@ +using System; +using Microsoft.VisualStudio.TestTools.UnitTesting; +using Numerics.Distributions; + +namespace Distributions.Univariate +{ + /// + /// Regression tests for the univariate distribution factory. + /// + [TestClass] + public class Test_UnivariateDistributionFactory + { + /// + /// Verifies that every defined enum member is either constructed exactly or rejected explicitly. + /// + [TestMethod] + public void EveryDefinedDistributionTypeIsHandledExplicitly() + { + foreach (UnivariateDistributionType type in Enum.GetValues(typeof(UnivariateDistributionType))) + { + if (type == UnivariateDistributionType.CompetingRisks || + type == UnivariateDistributionType.Mixture || + type == UnivariateDistributionType.UserDefined) + { + AssertThrows( + () => UnivariateDistributionFactory.CreateDistribution(type)); + } + else + { + var distribution = UnivariateDistributionFactory.CreateDistribution(type); + Assert.AreEqual(type, distribution.Type); + } + } + } + + /// + /// Verifies that undefined enum values fail fast rather than silently creating a deterministic distribution. + /// + [TestMethod] + public void UndefinedDistributionTypeThrowsArgumentOutOfRange() + { + AssertThrows( + () => UnivariateDistributionFactory.CreateDistribution((UnivariateDistributionType)int.MaxValue)); + } + + /// + /// Verifies that an action throws the requested exception type on every target framework. + /// + /// The exception type expected from . + /// The action expected to throw. + private static void AssertThrows(Action action) + where TException : Exception + { + try + { + action(); + } + catch (TException) + { + return; + } + catch (Exception exception) + { + Assert.Fail("Expected " + typeof(TException).Name + " but received " + exception.GetType().Name + "."); + return; + } + + Assert.Fail("Expected " + typeof(TException).Name + "."); + } + } +} diff --git a/Test_Numerics/Distributions/Univariate/Test_Weibull.cs b/Test_Numerics/Distributions/Univariate/Test_Weibull.cs index 3b3c2d04..b37d6354 100644 --- a/Test_Numerics/Distributions/Univariate/Test_Weibull.cs +++ b/Test_Numerics/Distributions/Univariate/Test_Weibull.cs @@ -1,4 +1,4 @@ -using System; +using System; using Microsoft.VisualStudio.TestTools.UnitTesting; using Numerics.Data.Statistics; using Numerics.Distributions; @@ -14,7 +14,7 @@ namespace Distributions.Univariate /// /// Haden Smith, USACE Risk Management Center, cole.h.smith@usace.army.mil /// Tiki Gonzalez, USACE Risk Management Center, julian.t.gonzalez@usace.army.mil - /// + /// /// /// /// References: @@ -50,7 +50,7 @@ public void Test_Weibull_MLE_Fit() /// /// /// - /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. + /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. /// /// /// Example 7.7.2 page 237. @@ -73,7 +73,7 @@ public void Test_Weibull_Quantile() /// /// /// - /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. + /// Reference: "Flood Frequency Analysis", A.R. Rao & K.H. Hamed, CRC Press, 2000. /// /// /// Example 7.2.3 page 240. diff --git a/Test_Numerics/Machine Learning/Supervised/Test_DecisionTree.cs b/Test_Numerics/Machine Learning/Supervised/Test_DecisionTree.cs index ef6dd210..fb087f34 100644 --- a/Test_Numerics/Machine Learning/Supervised/Test_DecisionTree.cs +++ b/Test_Numerics/Machine Learning/Supervised/Test_DecisionTree.cs @@ -1,4 +1,4 @@ -using Microsoft.VisualStudio.TestTools.UnitTesting; +using Microsoft.VisualStudio.TestTools.UnitTesting; using Numerics.Data; using Numerics.MachineLearning; using Numerics.Mathematics.LinearAlgebra; @@ -14,7 +14,7 @@ namespace MachineLearning /// /// /// Authors: - /// + /// /// Haden Smith, USACE Risk Management Center, cole.h.smith@usace.army.mil /// Tiki Gonzalez, USACE Risk Management Center, julian.t.gonzalez@usace.army.mil /// @@ -24,7 +24,7 @@ namespace MachineLearning public class Test_DecisionTree { /// - /// Decision Tree Classification is tested against the Iris dataset. + /// Decision Tree Classification is tested against the Iris dataset. /// [TestMethod] public void Test_DecisionTree_Iris() @@ -65,11 +65,11 @@ public void Test_DecisionTree_Iris() } /// - /// Testing the Decision Tree (DT) regression method. This test is mainly meant for demonstration. - /// I compare the DT regression against linear regression, and show the DT has worse performance against the observed data. + /// Testing the Decision Tree (DT) regression method. This test is mainly meant for demonstration. + /// I compare the DT regression against linear regression, and show the DT has worse performance against the observed data. /// Use a Random Forest to get better performance. /// - /// + /// /// /// [TestMethod] diff --git a/Test_Numerics/Machine Learning/Supervised/Test_GeneralizedLinearModel.cs b/Test_Numerics/Machine Learning/Supervised/Test_GeneralizedLinearModel.cs index 3adf17f4..d10f075b 100644 --- a/Test_Numerics/Machine Learning/Supervised/Test_GeneralizedLinearModel.cs +++ b/Test_Numerics/Machine Learning/Supervised/Test_GeneralizedLinearModel.cs @@ -1,4 +1,4 @@ -using Microsoft.VisualStudio.TestTools.UnitTesting; +using Microsoft.VisualStudio.TestTools.UnitTesting; using Numerics.Data; using Numerics.Functions; using Numerics.MachineLearning; @@ -18,7 +18,7 @@ namespace MachineLearning /// /// /// Authors: - /// + /// /// Haden Smith, USACE Risk Management Center, cole.h.smith@usace.army.mil /// /// @@ -28,7 +28,7 @@ public class Test_GeneralizedLinearModel { /// - /// Test simple linear regression. + /// Test simple linear regression. /// [TestMethod] public void Test_SimpleLinearRegression() @@ -51,7 +51,7 @@ public void Test_SimpleLinearRegression() var true_se = 0.6026; var true_sigA = 0.05569; var true_sigB = 0.04744; - // var true_r2 = 0.159; Never used below. + // var true_r2 = 0.159; Never used below. var true_df = 185; Assert.AreEqual(a, true_a, 1E-3); @@ -61,19 +61,19 @@ public void Test_SimpleLinearRegression() Assert.AreEqual(se, true_se, 1E-3); Assert.AreEqual(df, true_df); - // Test summary output table + // Test summary output table var summary = GLM.Summary(); for (int i = 0; i < summary.Count; i++) { //Debug.WriteLine(summary[i]); } - /** + /* Generalized Linear Model (187 obs, 2 parameters): Model for predicting Consumption: Parameters: Estimate Std. Error z value Pr(>|z|) - Intercept 0.54510 0.05569 9.789 < 1E-15 *** + Intercept 0.54510 0.05569 9.789 < 1E-15 *** β1 0.28057 0.04744 5.914 3.34E-009 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 @@ -90,7 +90,7 @@ Min 1Q Median 3Q Max /// - /// Test multiple linear regression. + /// Test multiple linear regression. /// [TestMethod] public void Test_MultipleLinearRegression() @@ -127,7 +127,7 @@ public void Test_MultipleLinearRegression() Assert.AreEqual(df, true_df); - // Test summary output table + // Test summary output table var summary = GLM.Summary(); for (int i = 0; i < summary.Count; i++) { @@ -135,15 +135,15 @@ public void Test_MultipleLinearRegression() //Debug.WriteLine(summary[i]); } - /** + /* Generalized Linear Model (187 obs, 5 parameters): Model for predicting Consumption: Parameters: Estimate Std. Error z value Pr(>|z|) Intercept 0.26729 0.03721 7.184 6.79E-013 *** - β1 0.71448 0.04219 16.934 < 1E-15 *** + β1 0.71448 0.04219 16.934 < 1E-15 *** β2 0.04588 0.02588 1.773 0.0762 . - β3 -0.04527 0.00278 -16.286 < 1E-15 *** + β3 -0.04527 0.00278 -16.286 < 1E-15 *** β4 -0.20484 0.10550 -1.942 0.0522 . --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 @@ -192,7 +192,7 @@ public void Test_Log() Assert.AreEqual(aic, true_aic, 1E-2); Assert.AreEqual(df, true_df); - // Test summary output table + // Test summary output table var summary = GLM.Summary(); for (int i = 0; i < summary.Count; i++) { @@ -200,14 +200,14 @@ public void Test_Log() Debug.WriteLine(summary[i]); } - /** + /* Generalized Linear Model (26 obs, 3 parameters): Model for predicting admits: Parameters: Estimate Std. Error z value Pr(>|z|) - Intercept 6.32270 0.00033 19,261.610 < 1E-15 *** - β1 0.00240 0.00000 5,917.117 < 1E-15 *** - β2 -0.00015 0.00000 -200.082 < 1E-15 *** + Intercept 6.32270 0.00033 19,261.610 < 1E-15 *** + β1 0.00240 0.00000 5,917.117 < 1E-15 *** + β2 -0.00015 0.00000 -200.082 < 1E-15 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 @@ -256,7 +256,7 @@ public void Test_Logistic() Assert.AreEqual(aic, true_aic, 1E-2); Assert.AreEqual(df, true_df); - // Test summary output table + // Test summary output table var summary = GLM.Summary(); for (int i = 0; i < summary.Count; i++) { @@ -264,7 +264,7 @@ public void Test_Logistic() //Debug.WriteLine(summary[i]); } - /** + /* Generalized Linear Model (400 obs, 4 parameters): Model for predicting admits: Parameters: @@ -321,7 +321,7 @@ public void Test_Probit() Assert.AreEqual(aic, true_aic, 1E-2); Assert.AreEqual(df, true_df); - // Test summary output table + // Test summary output table var summary = GLM.Summary(); for (int i = 0; i < summary.Count; i++) { @@ -329,7 +329,7 @@ public void Test_Probit() //Debug.WriteLine(summary[i]); } - /** + /* Generalized Linear Model (400 obs, 4 parameters): Model for predicting admits: Parameters: @@ -388,7 +388,7 @@ public void Test_LogLog() Assert.AreEqual(aic, true_aic, 1E-2); Assert.AreEqual(df, true_df); - // Test summary output table + // Test summary output table var summary = GLM.Summary(); for (int i = 0; i < summary.Count; i++) { @@ -396,7 +396,7 @@ public void Test_LogLog() Debug.WriteLine(summary[i]); } - /** + /* Generalized Linear Model (400 obs, 4 parameters): Model for predicting admits: Parameters: diff --git a/Test_Numerics/Machine Learning/Supervised/Test_RandomForest.cs b/Test_Numerics/Machine Learning/Supervised/Test_RandomForest.cs index fbd24229..491cfeaa 100644 --- a/Test_Numerics/Machine Learning/Supervised/Test_RandomForest.cs +++ b/Test_Numerics/Machine Learning/Supervised/Test_RandomForest.cs @@ -1,4 +1,4 @@ -using Microsoft.VisualStudio.TestTools.UnitTesting; +using Microsoft.VisualStudio.TestTools.UnitTesting; using Numerics.Data.Statistics; using Numerics.MachineLearning; using Numerics.Mathematics.LinearAlgebra; @@ -14,7 +14,7 @@ namespace MachineLearning /// /// /// Authors: - /// + /// /// Haden Smith, USACE Risk Management Center, cole.h.smith@usace.army.mil /// /// @@ -23,7 +23,7 @@ namespace MachineLearning public class Test_RandomForest { /// - /// Random Forest Classification is tested against the Iris dataset. + /// Random Forest Classification is tested against the Iris dataset. /// [TestMethod] public void Test_RandomForest_Iris() @@ -64,10 +64,10 @@ public void Test_RandomForest_Iris() } /// - /// Testing the Random Forest (RF) regression method. This test is mainly meant for demonstration. - /// I compare the RF regression against linear regression, and show the RF has better performance against the observed data. + /// Testing the Random Forest (RF) regression method. This test is mainly meant for demonstration. + /// I compare the RF regression against linear regression, and show the RF has better performance against the observed data. /// - /// + /// /// /// [TestMethod] diff --git a/Test_Numerics/Machine Learning/Supervised/Test_kNN.cs b/Test_Numerics/Machine Learning/Supervised/Test_kNN.cs index 032e244e..6dc8f883 100644 --- a/Test_Numerics/Machine Learning/Supervised/Test_kNN.cs +++ b/Test_Numerics/Machine Learning/Supervised/Test_kNN.cs @@ -1,4 +1,4 @@ -using Microsoft.VisualStudio.TestTools.UnitTesting; +using Microsoft.VisualStudio.TestTools.UnitTesting; using Numerics.MachineLearning; using Numerics.Mathematics.LinearAlgebra; using System.Collections.Generic; @@ -14,7 +14,7 @@ namespace MachineLearning /// /// /// Authors: - /// + /// /// Haden Smith, USACE Risk Management Center, cole.h.smith@usace.army.mil /// Tiki Gonzalez, USACE Risk Management Center, julian.t.gonzalez@usace.army.mil /// @@ -67,7 +67,7 @@ public void Test_kNN_Iris() /// /// Testing classification method in kNN. Expected a classification of 0. /// - /// + /// /// /// [TestMethod()] @@ -88,10 +88,10 @@ public void Test_kNN_Classification() } /// - /// Testing regression method in kNN. This test is mainly meant for demonstration. + /// Testing regression method in kNN. This test is mainly meant for demonstration. /// I compare kNN regression against linear regression, and show kNN has better performance against the observed data. /// - /// + /// /// /// [TestMethod] @@ -130,7 +130,7 @@ public void Test_kNN_Regression() // Get R-Squared of predictions var knnR2 = GoodnessOfFit.RSquared(Y_test.Array, knnPredict); var lmR2 = GoodnessOfFit.RSquared(Y_test.Array, lmPredict); - + // kNN is better Assert.IsGreaterThan(lmR2, knnR2 ); diff --git a/Test_Numerics/Mathematics/Integration/Integrands.cs b/Test_Numerics/Mathematics/Integration/Integrands.cs index 46487ffe..addec3c2 100644 --- a/Test_Numerics/Mathematics/Integration/Integrands.cs +++ b/Test_Numerics/Mathematics/Integration/Integrands.cs @@ -55,7 +55,8 @@ public static double FXXX(double x) /// /// Test function. The integral of Pi. Should equal ~3.14 /// - /// Array of values. + /// The x-coordinate. + /// The y-coordinate. public static double PI2D(double x, double y) { return (x * x + y * y < 1) ? 1 : 0; @@ -81,7 +82,14 @@ public static double GSL(double[] x) return A / (1.0 - Math.Cos(x[0]) * Math.Cos(x[1]) * Math.Cos(x[2])); } + /// + /// Mean values for the 20-dimensional sum-of-normals integration fixture. + /// public static double[] mu20 = new double[] { 10, 30, 17, 99, 68, 26, 35, 55, 13, 59, 12, 28, 49, 54, 20, 47, 12, 76, 70, 57 }; + + /// + /// Standard deviation values for the 20-dimensional sum-of-normals integration fixture. + /// public static double[] sigma20 = new double[] { 2, 15, 5, 14, 7, 24, 29, 22, 22, 1, 3, 28, 19, 18, 4, 24, 23, 26, 26, 19 }; /// diff --git a/Test_Numerics/Mathematics/Integration/Test_Integration.cs b/Test_Numerics/Mathematics/Integration/Test_Integration.cs index 0e4a7966..926e6750 100644 --- a/Test_Numerics/Mathematics/Integration/Test_Integration.cs +++ b/Test_Numerics/Mathematics/Integration/Test_Integration.cs @@ -1,4 +1,4 @@ -using Microsoft.VisualStudio.TestTools.UnitTesting; +using Microsoft.VisualStudio.TestTools.UnitTesting; namespace Mathematics.Integration { @@ -56,7 +56,7 @@ public void Test_TrapezoidalRule() /// - /// Test Simpsons Rule Method. + /// Test Simpsons Rule Method. /// [TestMethod()] public void Test_SimpsonsRule() @@ -68,7 +68,7 @@ public void Test_SimpsonsRule() } /// - /// Midpoint Method. + /// Midpoint Method. /// [TestMethod()] public void Test_MidPoint() @@ -79,4 +79,4 @@ public void Test_MidPoint() } } -} \ No newline at end of file +} diff --git a/Test_Numerics/Mathematics/Linear Algebra/Test_CholeskyDecomposition.cs b/Test_Numerics/Mathematics/Linear Algebra/Test_CholeskyDecomposition.cs index 25d5722b..25386978 100644 --- a/Test_Numerics/Mathematics/Linear Algebra/Test_CholeskyDecomposition.cs +++ b/Test_Numerics/Mathematics/Linear Algebra/Test_CholeskyDecomposition.cs @@ -1,4 +1,4 @@ -using System; +using System; using Microsoft.VisualStudio.TestTools.UnitTesting; using Numerics.Mathematics.LinearAlgebra; @@ -10,7 +10,7 @@ namespace Mathematics.LinearAlgebra /// /// /// Authors: - /// + /// /// Haden Smith, USACE Risk Management Center, cole.h.smith@usace.army.mil /// Tiki Gonzalez, USACE Risk Management Center, julian.t.gonzalez@usace.army.mil /// @@ -62,14 +62,13 @@ public void Test_CholeskyDecomp() /// /// Test Cholesky Decomposition method with known values from textbook. - /// /// - - + + [TestMethod()] public void Test_CholeskyDecompVals() { - + var A = new Matrix(3); A[0, 0] = 16d; A[0, 1] = 4d; @@ -94,7 +93,7 @@ public void Test_CholeskyDecompVals() L[2, 1] = -3d; L[2, 2] = 3d; - //Test Decomposition + //Test Decomposition var prod = L; for (int i = 0; i < prod.NumberOfRows; i++) { @@ -121,7 +120,7 @@ public void Test_LogDeterminant() var chol = new CholeskyDecomposition(A); // Test Determinant - double true_det = 6.356108; //Math.Log(576) + double true_det = 6.356108; //Math.Log(576) double det = chol.LogDeterminant(); Assert.AreEqual(det, true_det, 0.0001d); } @@ -182,7 +181,7 @@ public void Test_Solve() A[2, 2] = 22d; var chol = new CholeskyDecomposition(A); - //Test Solve + //Test Solve var B = new Vector(new[] { 16d, 18d, -22d }); var trueX = new Vector(new[] { 1d, 2d, -1d }); var X = chol.Solve(B); @@ -192,7 +191,7 @@ public void Test_Solve() /// /// Test Forward-substitution method - /// + /// [TestMethod()] public void Test_Forward() { @@ -207,7 +206,7 @@ public void Test_Forward() A[2, 1] = -4d; A[2, 2] = 22d; var chol = new CholeskyDecomposition(A); - + var b = new Vector(new double[] { 16d, 18d, -22d }); var true_y = new double[] { 4d, 7d, -3d }; var y = chol.Forward(b); @@ -217,7 +216,7 @@ public void Test_Forward() /// /// Tests Back-substitution method - /// + /// [TestMethod()] public void Test_Back() { @@ -232,14 +231,14 @@ public void Test_Back() A[2, 1] = -4d; A[2, 2] = 22d; var chol = new CholeskyDecomposition(A); - + var Y = new Vector(new double[] { 4d, 7d, -3d }); var right_x = new double[] { 1d, 2d, -1d }; var x = chol.Backward(Y); for (int i = 0; i < x.Length; i++) Assert.AreEqual(x[i], right_x[i], 0.0001d); } - + } } diff --git a/Test_Numerics/Mathematics/Linear Algebra/Test_GaussJordanElimination.cs b/Test_Numerics/Mathematics/Linear Algebra/Test_GaussJordanElimination.cs index 4e3d166a..d773749e 100644 --- a/Test_Numerics/Mathematics/Linear Algebra/Test_GaussJordanElimination.cs +++ b/Test_Numerics/Mathematics/Linear Algebra/Test_GaussJordanElimination.cs @@ -1,4 +1,4 @@ -using Microsoft.VisualStudio.TestTools.UnitTesting; +using Microsoft.VisualStudio.TestTools.UnitTesting; using Numerics.Mathematics.LinearAlgebra; namespace Mathematics.LinearAlgebra @@ -9,7 +9,7 @@ namespace Mathematics.LinearAlgebra /// /// /// Authors: - /// + /// /// Haden Smith, USACE Risk Management Center, cole.h.smith@usace.army.mil /// Tiki Gonzalez, USACE Risk Management Center, julian.t.gonzalez@usace.army.mil /// @@ -18,6 +18,9 @@ namespace Mathematics.LinearAlgebra [TestClass] public class Test_GaussJordanElimination { + /// + /// Tests Gauss-Jordan elimination against a known solution. + /// [TestMethod()] public void Test_GaussJordanElim() { @@ -32,8 +35,8 @@ public void Test_GaussJordanElim() Assert.AreEqual(true_IA[i, j], A[i, j]); } - /// Recreated Gauss Jordan test in R to compare the inverted A matrices. - /// I utilized library(matlib), gaussianElimination(), and inv() functions. Test passed. + // Recreated Gauss Jordan test in R to compare the inverted A matrices. + // I utilized library(matlib), gaussianElimination(), and inv() functions. Test passed. } } diff --git a/Test_Numerics/Mathematics/Linear Algebra/Test_LUDecomposition.cs b/Test_Numerics/Mathematics/Linear Algebra/Test_LUDecomposition.cs index 0e95f9b8..e282bae5 100644 --- a/Test_Numerics/Mathematics/Linear Algebra/Test_LUDecomposition.cs +++ b/Test_Numerics/Mathematics/Linear Algebra/Test_LUDecomposition.cs @@ -1,4 +1,4 @@ -using Microsoft.VisualStudio.TestTools.UnitTesting; +using Microsoft.VisualStudio.TestTools.UnitTesting; using Numerics.Mathematics.LinearAlgebra; namespace Mathematics.LinearAlgebra @@ -9,7 +9,7 @@ namespace Mathematics.LinearAlgebra /// /// /// Authors: - /// + /// /// Haden Smith, USACE Risk Management Center, cole.h.smith@usace.army.mil /// Tiki Gonzalez, USACE Risk Management Center, julian.t.gonzalez@usace.army.mil /// @@ -67,11 +67,11 @@ public void Test_LUDecomp() } } - /// Replicated LU decomposition in R and compared with decomposition results. - /// Used lu() function for recreation. Test passed. + // Replicated LU decomposition in R and compared with decomposition results. + // Used lu() function for recreation. Test passed. } /// - /// Testing Solve with vector input + /// Testing Solve with vector input /// [TestMethod()] public void Test_SolveVector() { @@ -89,7 +89,7 @@ public void Test_SolveVector() { var B = new Vector(new[] { 6d, -4, 27d }); var lu = new LUDecomposition(A); - + var true_x = new[] { 5d, 3d, -2 }; var x = lu.Solve(B); for (int i = 0; i < x.Length; i++) @@ -97,7 +97,7 @@ public void Test_SolveVector() { } /// - /// Testing Solve with matrix input + /// Testing Solve with matrix input /// [TestMethod()] public void Test_SolveMatrix() { @@ -130,9 +130,9 @@ public void Test_SolveMatrix() { for (int i = 0; i < matX.NumberOfRows; i++) { Assert.AreEqual(matX[i,0], true_matX[i,0], 0.0001d); - + } - + } /// diff --git a/Test_Numerics/Mathematics/Linear Algebra/Test_Matrix.cs b/Test_Numerics/Mathematics/Linear Algebra/Test_Matrix.cs index 7dc9be74..ed221acc 100644 --- a/Test_Numerics/Mathematics/Linear Algebra/Test_Matrix.cs +++ b/Test_Numerics/Mathematics/Linear Algebra/Test_Matrix.cs @@ -1,15 +1,15 @@ -using Microsoft.VisualStudio.TestTools.UnitTesting; +using Microsoft.VisualStudio.TestTools.UnitTesting; using Numerics.Mathematics.LinearAlgebra; namespace Mathematics.LinearAlgebra -{ +{ /// /// A class unit testing all Matrix Operations. /// /// /// /// Authors: - /// + /// /// Haden Smith, USACE Risk Management Center, cole.h.smith@usace.army.mil /// Tiki Gonzalez, USACE Risk Management Center, julian.t.gonzalez@usace.army.mil /// @@ -35,33 +35,33 @@ public void Test_StructuralMethods() A[2, 1] = 8d; A[2, 2] = 9d; - /// - /// Testing .Row - /// + // + // Testing .Row + // var true_row = new[] { 4d, 5d, 6d }; var row = A.Row(1); for (int i = 0; i < row.Length; i++) Assert.AreEqual(row[i], true_row[i]); - /// - /// Testing .Column - /// + // + // Testing .Column + // var true_col = new[] { 2d, 5d, 8d }; var col = A.Column(1); for (int i = 0; i < col.Length; i++) Assert.AreEqual(col[i], true_col[i]); - /// - ///Testing .Diagonal - /// + // + //Testing .Diagonal + // var true_diag = new[] { 1d, 5d, 9d }; var diag = A.Diagonal(); for (int i = 0; i < diag.Length; i++) Assert.AreEqual(diag[i], true_diag[i]); - /// - ///Testing upper triangular - /// + // + //Testing upper triangular + // var true_upT = new Matrix(3); true_upT[0, 0] = 1d; true_upT[0, 1] = 2d; @@ -76,9 +76,9 @@ public void Test_StructuralMethods() Assert.AreEqual(upT[i, j], true_upT[i, j]); } - /// - ///Testing lower triangular - /// + // + //Testing lower triangular + // var true_lowT = new Matrix(3); true_lowT[0, 0] = 1d; true_lowT[1, 0] = 4d; @@ -330,7 +330,7 @@ public void Test_Diagonal() Assert.AreEqual(result[i, j], true_result[i, j]); } } - + /// /// Sum all elements in the matrix /// @@ -479,7 +479,7 @@ public void Test_MultiplyByVector() } /// - /// Multiply matrix by scalar. + /// Multiply matrix by scalar. /// [TestMethod()] public void Test_MultiplyByScalar() diff --git a/Test_Numerics/Mathematics/Linear Algebra/Test_SingularValueDecomp.cs b/Test_Numerics/Mathematics/Linear Algebra/Test_SingularValueDecomp.cs index 43eccc36..46d675e7 100644 --- a/Test_Numerics/Mathematics/Linear Algebra/Test_SingularValueDecomp.cs +++ b/Test_Numerics/Mathematics/Linear Algebra/Test_SingularValueDecomp.cs @@ -1,4 +1,4 @@ -using System; +using System; using Microsoft.VisualStudio.TestTools.UnitTesting; using Numerics.Mathematics.LinearAlgebra; @@ -10,7 +10,7 @@ namespace Mathematics.LinearAlgebra /// /// /// Authors: - /// + /// /// Haden Smith, USACE Risk Management Center, cole.h.smith@usace.army.mil /// Tiki Gonzalez, USACE Risk Management Center, julian.t.gonzalez@usace.army.mil /// @@ -64,7 +64,7 @@ public void Test_Rank() } /// - /// Tests nullity function. + /// Tests nullity function. /// [TestMethod()] public void Test_Nullity() @@ -86,7 +86,7 @@ public void Test_Nullity() } /// - /// Tests known Range of matrix. + /// Tests known Range of matrix. /// The function should return the column space of the matrix. /// [TestMethod()] @@ -129,7 +129,7 @@ public void Test_Range() } /// - /// Tests known Nullspace of a matrix + /// Tests known Nullspace of a matrix /// [TestMethod()] public void Test_Nullspace() @@ -175,7 +175,7 @@ public void Test_SolveVector() var B = new Vector(new[] { 7d, -2d, 0 }); - //Testing Solve with vector input + //Testing Solve with vector input var true_x = new[] { 1d, -1d, 2 }; var x = svd.Solve(B); for (int i = 0; i < x.Length; i++) @@ -183,7 +183,7 @@ public void Test_SolveVector() } /// - /// Evaluating A*x = B with a B matrix input + /// Evaluating A*x = B with a B matrix input /// [TestMethod()] public void Test_SolveMatrix() @@ -265,13 +265,13 @@ public void Test_Decompose() Assert.AreEqual(V[i, j], svd.V[i, j], 0.0001d); } - ///Tested with the Test_Decompose() matrix to see if U and V were equal - ///svd() function assigns switched left and right singular vectors, however - ///the overall decomposition is the same for both functions. Test passed. + //Tested with the Test_Decompose() matrix to see if U and V were equal + //svd() function assigns switched left and right singular vectors, however + //the overall decomposition is the same for both functions. Test passed. } /// - /// Testing log determinant of decomposed matrix + /// Testing log determinant of decomposed matrix /// [TestMethod()] public void Test_LogDeterminant() @@ -289,7 +289,7 @@ public void Test_LogDeterminant() var svd = new SingularValueDecomposition(A); // Test Determinant - double true_det = 6.356108; //Math.Log(576) + double true_det = 6.356108; //Math.Log(576) double det = svd.LogDeterminant(); Assert.AreEqual(det, true_det, 0.0001d); } @@ -313,7 +313,7 @@ public void Test_LogPseudoDeterminant() var svd = new SingularValueDecomposition(A); // Test Determinant - double true_det = 6.356108; //Math.Log(576) + double true_det = 6.356108; //Math.Log(576) double det = svd.LogPseudoDeterminant(); Assert.AreEqual(det, true_det, 0.0001d); } diff --git a/Test_Numerics/Mathematics/ODE Solvers/Test_RungeKutta.cs b/Test_Numerics/Mathematics/ODE Solvers/Test_RungeKutta.cs index bc652100..241d037a 100644 --- a/Test_Numerics/Mathematics/ODE Solvers/Test_RungeKutta.cs +++ b/Test_Numerics/Mathematics/ODE Solvers/Test_RungeKutta.cs @@ -79,6 +79,7 @@ public void Test_4thRK() /// /// Test the second fourth order Runge-Kutta method that allows you to specify the time step size (dt) instead of the number of time steps + /// [TestMethod] public void Test_4thRK_2() { diff --git a/Test_Numerics/Mathematics/Optimization/Dynamic/BinaryHeapTesting.cs b/Test_Numerics/Mathematics/Optimization/Dynamic/BinaryHeapTesting.cs index 46e5ad8e..b7028b33 100644 --- a/Test_Numerics/Mathematics/Optimization/Dynamic/BinaryHeapTesting.cs +++ b/Test_Numerics/Mathematics/Optimization/Dynamic/BinaryHeapTesting.cs @@ -1,4 +1,4 @@ -using System; +using System; using System.Collections.Generic; using System.Diagnostics; using Microsoft.VisualStudio.TestTools.UnitTesting; @@ -7,6 +7,9 @@ namespace Mathematics.Optimization { + /// + /// Tests binary heap behavior used by dynamic optimization routines. + /// [TestClass] public class BinaryHeapTesting { @@ -119,7 +122,7 @@ public void HeapTest4() } /// - /// Decrease key for ordering + /// Decrease key for ordering /// [TestMethod] public void DecreaseKeyTest() @@ -196,4 +199,4 @@ public void OrderingWithNegativeWeightsTest() } -} \ No newline at end of file +} diff --git a/Test_Numerics/Mathematics/Optimization/Dynamic/DijkstraTesting.cs b/Test_Numerics/Mathematics/Optimization/Dynamic/DijkstraTesting.cs index 2185a832..030c673b 100644 --- a/Test_Numerics/Mathematics/Optimization/Dynamic/DijkstraTesting.cs +++ b/Test_Numerics/Mathematics/Optimization/Dynamic/DijkstraTesting.cs @@ -1,8 +1,11 @@ - + using Numerics.Mathematics.Optimization; namespace Mathematics.Optimization { + /// + /// Tests shortest-path routing behavior for Dijkstra networks. + /// [TestClass] public class ShortestPathTesting { @@ -32,7 +35,7 @@ public void SimpleEdgeGraphCost() } /// - /// Simple network run, testing to see if algorithm chooses + /// Simple network run, testing to see if algorithm chooses /// the lowest cost path as it should. /// [TestMethod] @@ -223,7 +226,7 @@ public void MultipleDestSharedPath() [TestMethod] public void DisconnectedComponent() { - // Graph: + // Graph: // 0 - 1 2 - 3 var edges = new List { diff --git a/Test_Numerics/Mathematics/Optimization/Global/Test_MLSL.cs b/Test_Numerics/Mathematics/Optimization/Global/Test_MLSL.cs index 6f21fadb..4fcd361e 100644 --- a/Test_Numerics/Mathematics/Optimization/Global/Test_MLSL.cs +++ b/Test_Numerics/Mathematics/Optimization/Global/Test_MLSL.cs @@ -1,4 +1,4 @@ -using Microsoft.VisualStudio.TestTools.UnitTesting; +using Microsoft.VisualStudio.TestTools.UnitTesting; using Numerics.Mathematics.Optimization; namespace Mathematics.Optimization @@ -9,7 +9,7 @@ namespace Mathematics.Optimization /// /// /// Authors: - /// + /// /// Haden Smith, USACE Risk Management Center, cole.h.smith@usace.army.mil /// /// @@ -174,7 +174,7 @@ public void Test_McCormick() /// /// Test the MLSL algorithm with the Beale Function - /// [TestMethod] public void Test_Beale() { diff --git a/Test_Numerics/Mathematics/Optimization/Global/Test_MultiStart.cs b/Test_Numerics/Mathematics/Optimization/Global/Test_MultiStart.cs index 78509209..f02bd6f4 100644 --- a/Test_Numerics/Mathematics/Optimization/Global/Test_MultiStart.cs +++ b/Test_Numerics/Mathematics/Optimization/Global/Test_MultiStart.cs @@ -1,4 +1,4 @@ -using Microsoft.VisualStudio.TestTools.UnitTesting; +using Microsoft.VisualStudio.TestTools.UnitTesting; using Numerics.Mathematics.Optimization; namespace Mathematics.Optimization @@ -9,7 +9,7 @@ namespace Mathematics.Optimization /// /// /// Authors: - /// + /// /// Haden Smith, USACE Risk Management Center, cole.h.smith@usace.army.mil /// /// @@ -174,7 +174,7 @@ public void Test_McCormick() /// /// Test the MS algorithm with the Beale Function - /// [TestMethod] public void Test_Beale() { diff --git a/Test_Numerics/Mathematics/Optimization/Local/Test_BrentSearch.cs b/Test_Numerics/Mathematics/Optimization/Local/Test_BrentSearch.cs index 45cd0b86..b795d9b1 100644 --- a/Test_Numerics/Mathematics/Optimization/Local/Test_BrentSearch.cs +++ b/Test_Numerics/Mathematics/Optimization/Local/Test_BrentSearch.cs @@ -3,6 +3,9 @@ namespace Mathematics.Optimization { + using System; + using System.Collections.Generic; + /// /// Unit tests for the Brent optimization algorithm /// @@ -63,5 +66,199 @@ public void Test_DeJong() Assert.AreEqual(X, trueX, 1E-4); } + /// + /// Verifies geometric expansion to the right and the exact trial sequence. + /// + [TestMethod] + public void Test_Bracket_GeometricExpansionRight() + { + var evaluations = new List(); + var solver = new BrentSearch(x => + { + evaluations.Add(x); + return (x - 10d) * (x - 10d); + }, 0d, 1d); + + solver.Bracket(1d, 2d); + + CollectionAssert.AreEqual(new[] { 0d, 1d, 2d, 4d, 8d, 16d }, evaluations.ToArray()); + Assert.AreEqual(4d, solver.LowerBound, 0d); + Assert.AreEqual(16d, solver.UpperBound, 0d); + } + + /// + /// Verifies that a positive initial step reverses direction and expands geometrically when the minimum is to the left. + /// + [TestMethod] + public void Test_Bracket_GeometricExpansionLeft() + { + var evaluations = new List(); + var solver = new BrentSearch(x => + { + evaluations.Add(x); + return (x + 10d) * (x + 10d); + }, 0d, 1d); + + solver.Bracket(1d, 2d); + + CollectionAssert.AreEqual(new[] { 0d, 1d, -1d, -3d, -7d, -15d }, evaluations.ToArray()); + Assert.AreEqual(-15d, solver.LowerBound, 0d); + Assert.AreEqual(-3d, solver.UpperBound, 0d); + } + + /// + /// Verifies that a caller-provided expansion factor controls the geometric trial sequence. + /// + [TestMethod] + public void Test_Bracket_UsesCustomExpansionFactor() + { + var evaluations = new List(); + var solver = new BrentSearch(x => + { + evaluations.Add(x); + return (x - 10d) * (x - 10d); + }, 0d, 1d); + + solver.Bracket(1d, 3d); + + CollectionAssert.AreEqual(new[] { 0d, 1d, 2d, 5d, 14d, 41d }, evaluations.ToArray()); + Assert.AreEqual(5d, solver.LowerBound, 0d); + Assert.AreEqual(41d, solver.UpperBound, 0d); + } + + /// + /// Verifies that geometric bracketing finds a distant minimum with logarithmic objective work and remains compatible with minimization. + /// + [TestMethod] + public void Test_Bracket_DistantMinimumReducesEvaluations() + { + int evaluations = 0; + var solver = new BrentSearch(x => + { + evaluations++; + return (x - 1000d) * (x - 1000d); + }, 0d, 1d); + + solver.Bracket(0.1d, 2d); + int bracketEvaluations = evaluations; + + Assert.IsLessThanOrEqualTo(20, bracketEvaluations, $"Expected at most 20 bracket evaluations, but observed {bracketEvaluations}."); + Assert.IsTrue(solver.LowerBound <= 1000d && solver.UpperBound >= 1000d); + solver.Minimize(); + Assert.AreEqual(1000d, solver.BestParameterSet.Values[0], 1E-4); + Assert.AreEqual(0d, solver.BestParameterSet.Fitness, 1E-8); + } + + /// + /// Verifies that a non-strict bracket terminates immediately for a flat objective. + /// + [TestMethod] + public void Test_Bracket_PlateauTerminates() + { + int evaluations = 0; + var solver = new BrentSearch(x => + { + evaluations++; + return 1d; + }, 0d, 1d); + + solver.Bracket(1d, 2d); + + Assert.AreEqual(3, evaluations); + Assert.AreEqual(0d, solver.LowerBound, 0d); + Assert.AreEqual(2d, solver.UpperBound, 0d); + } + + /// + /// Verifies bounded failure for a monotone objective and transactional preservation of the original bounds. + /// + [TestMethod] + public void Test_Bracket_MonotoneObjectiveReachesIterationLimit() + { + var solver = new BrentSearch(x => -x, 0d, 1d) + { + MaxIterations = 10, + ReportFailure = false + }; + + solver.Bracket(1d, 2d); + + Assert.AreEqual(OptimizationStatus.MaximumIterationsReached, solver.Status); + Assert.AreEqual(0d, solver.LowerBound, 0d); + Assert.AreEqual(1d, solver.UpperBound, 0d); + + var throwingSolver = new BrentSearch(x => -x, 0d, 1d) { MaxIterations = 10 }; + var exception = AssertThrows(() => throwingSolver.Bracket(1d, 2d)); + Assert.AreEqual(nameof(Optimizer.MaxIterations), exception.ParamName); + } + + /// + /// Verifies validation of the public bracketing step and expansion-factor contract. + /// + [TestMethod] + public void Test_Bracket_RejectsInvalidInputs() + { + foreach (double invalidStep in new[] { 0d, double.NaN, double.NegativeInfinity, double.PositiveInfinity }) + { + var solver = new BrentSearch(x => x * x, 0d, 1d); + var exception = AssertThrows(() => solver.Bracket(invalidStep, 2d)); + Assert.AreEqual("s", exception.ParamName); + } + + foreach (double invalidExpansion in new[] { -1d, 0d, 1d, double.NaN, double.NegativeInfinity, double.PositiveInfinity }) + { + var solver = new BrentSearch(x => x * x, 0d, 1d); + var exception = AssertThrows(() => solver.Bracket(1d, invalidExpansion)); + Assert.AreEqual("k", exception.ParamName); + } + } + + /// + /// Verifies deterministic failure for NaN objectives and coordinate overflow without changing the original bounds. + /// + [TestMethod] + public void Test_Bracket_RejectsNonFiniteSearchState() + { + var nanSolver = new BrentSearch(x => double.NaN, 0d, 1d) { ReportFailure = false }; + nanSolver.Bracket(); + Assert.AreEqual(OptimizationStatus.Failure, nanSolver.Status); + Assert.AreEqual(0d, nanSolver.LowerBound, 0d); + Assert.AreEqual(1d, nanSolver.UpperBound, 0d); + + var overflowSolver = new BrentSearch(x => -x, double.MaxValue, double.MaxValue) { ReportFailure = false }; + overflowSolver.Bracket(double.MaxValue, 2d); + Assert.AreEqual(OptimizationStatus.Failure, overflowSolver.Status); + Assert.AreEqual(double.MaxValue, overflowSolver.LowerBound); + Assert.AreEqual(double.MaxValue, overflowSolver.UpperBound); + + var throwingSolver = new BrentSearch(x => -x, double.MaxValue, double.MaxValue); + AssertThrows(() => throwingSolver.Bracket(double.MaxValue, 2d)); + } + + /// + /// Executes an action and returns the expected exception. + /// + /// The expected exception type. + /// The action expected to throw. + /// The exception thrown by . + private static TException AssertThrows(Action action) where TException : Exception + { + try + { + action(); + } + catch (TException exception) + { + return exception; + } + catch (Exception exception) + { + Assert.Fail($"Expected {typeof(TException).Name}, but observed {exception.GetType().Name}."); + } + + Assert.Fail($"Expected {typeof(TException).Name}, but no exception was thrown."); + throw new InvalidOperationException("The exception assertion did not terminate as expected."); + } + } } diff --git a/Test_Numerics/Mathematics/Optimization/Local/Test_Powell.cs b/Test_Numerics/Mathematics/Optimization/Local/Test_Powell.cs index b7bb5f0e..b598290f 100644 --- a/Test_Numerics/Mathematics/Optimization/Local/Test_Powell.cs +++ b/Test_Numerics/Mathematics/Optimization/Local/Test_Powell.cs @@ -187,5 +187,29 @@ public void Test_Beale() Assert.AreEqual(y, validY, 1E-4); } + /// + /// Verifies that Powell's Brent line search reaches a distant minimum without linear bracketing work. + /// + [TestMethod] + public void Test_DistantLineMinimumUsesGeometricBracket() + { + var solver = new Powell( + x => (x[0] - 50d) * (x[0] - 50d), + 1, + new[] { 0d }, + new[] { -100d }, + new[] { 100d }) + { + ComputeHessian = false, + RecordTraces = false + }; + + solver.Minimize(); + + Assert.AreEqual(50d, solver.BestParameterSet.Values[0], 1E-4); + Assert.AreEqual(0d, solver.BestParameterSet.Fitness, 1E-8); + Assert.IsLessThan(200, solver.FunctionEvaluations, $"Expected fewer than 200 objective evaluations, but observed {solver.FunctionEvaluations}."); + } + } } diff --git a/Test_Numerics/Mathematics/Root Finding/Test_Bisection.cs b/Test_Numerics/Mathematics/Root Finding/Test_Bisection.cs index 4000c839..deb14c87 100644 --- a/Test_Numerics/Mathematics/Root Finding/Test_Bisection.cs +++ b/Test_Numerics/Mathematics/Root Finding/Test_Bisection.cs @@ -1,4 +1,4 @@ -using Microsoft.VisualStudio.TestTools.UnitTesting; +using Microsoft.VisualStudio.TestTools.UnitTesting; using Numerics.Mathematics.RootFinding; using System; @@ -10,7 +10,7 @@ namespace Mathematics.RootFinding /// /// /// Authors: - /// + /// /// Haden Smith, USACE Risk Management Center, cole.h.smith@usace.army.mil /// Tiki Gonzalez, USACE Risk Management Center, julian.t.gonzalez@usace.army.mil /// @@ -90,6 +90,9 @@ public void Test_Power() Assert.AreEqual(X, trueX, 1E-5); } + /// + /// Tests the first bisection edge case. + /// [TestMethod()] public void Test_BisectionEdge1() { @@ -102,6 +105,9 @@ public void Test_BisectionEdge1() }); } + /// + /// Tests the second bisection edge case. + /// [TestMethod()] public void Test_BisectionEdge2() { @@ -114,6 +120,9 @@ public void Test_BisectionEdge2() }); } + /// + /// Tests the third bisection edge case. + /// [TestMethod()] public void Test_BisectionEdge3() { diff --git a/Test_Numerics/Mathematics/Root Finding/Test_Bracket.cs b/Test_Numerics/Mathematics/Root Finding/Test_Bracket.cs index b2522d4a..5d43390c 100644 --- a/Test_Numerics/Mathematics/Root Finding/Test_Bracket.cs +++ b/Test_Numerics/Mathematics/Root Finding/Test_Bracket.cs @@ -1,4 +1,4 @@ -using Microsoft.VisualStudio.TestTools.UnitTesting; +using Microsoft.VisualStudio.TestTools.UnitTesting; using Numerics.Mathematics.RootFinding; using System; @@ -10,7 +10,7 @@ namespace Mathematics.RootFinding /// /// /// Authors: - /// + /// /// Haden Smith, USACE Risk Management Center, cole.h.smith@usace.army.mil /// Tiki Gonzalez, USACE Risk Management Center, julian.t.gonzalez@usace.army.mil /// @@ -21,7 +21,7 @@ public class Test_Bracket { /// - /// Test with a quadratic function. + /// Test with a quadratic function. /// [TestMethod()] public void Test_Quadratic() @@ -33,7 +33,7 @@ public void Test_Quadratic() } /// - /// Test with a cubic function. + /// Test with a cubic function. /// [TestMethod()] public void Test_Cubic() @@ -45,7 +45,7 @@ public void Test_Cubic() } /// - /// Test with a trigonometric function. + /// Test with a trigonometric function. /// [TestMethod()] public void Test_Trigonometric() @@ -99,4 +99,4 @@ public void Test_BracketEdge() }); } } -} \ No newline at end of file +} diff --git a/Test_Numerics/Mathematics/Root Finding/Test_Secant.cs b/Test_Numerics/Mathematics/Root Finding/Test_Secant.cs index 761ea095..6cdbb60d 100644 --- a/Test_Numerics/Mathematics/Root Finding/Test_Secant.cs +++ b/Test_Numerics/Mathematics/Root Finding/Test_Secant.cs @@ -1,4 +1,4 @@ -using Microsoft.VisualStudio.TestTools.UnitTesting; +using Microsoft.VisualStudio.TestTools.UnitTesting; using Numerics.Mathematics.RootFinding; using System; @@ -10,7 +10,7 @@ namespace Mathematics.RootFinding /// /// /// Authors: - /// + /// /// Haden Smith, USACE Risk Management Center, cole.h.smith@usace.army.mil /// Tiki Gonzalez, USACE Risk Management Center, julian.t.gonzalez@usace.army.mil /// @@ -60,7 +60,7 @@ public void Test_Exponential() } /// - /// Testing edge case where the function is discontinuous within the interval + /// Testing edge case where the function is discontinuous within the interval /// public double Undefined(double x) { @@ -68,6 +68,9 @@ public double Undefined(double x) return F; } + /// + /// Tests secant method edge-case behavior. + /// [TestMethod()] public void Test_Edge() { diff --git a/Test_Numerics/Mathematics/Special Functions/Test_Bessel.cs b/Test_Numerics/Mathematics/Special Functions/Test_Bessel.cs index a11539e0..11696ca5 100644 --- a/Test_Numerics/Mathematics/Special Functions/Test_Bessel.cs +++ b/Test_Numerics/Mathematics/Special Functions/Test_Bessel.cs @@ -252,7 +252,7 @@ public void Test_Jn_ForwardRecurrence() } /// - /// Test Jn using the Miller's downward recurrence path (x <= n). + /// Test Jn using the Miller's downward recurrence path (x <= n). /// [TestMethod] public void Test_Jn_MillerRecurrence() diff --git a/Test_Numerics/Mathematics/Special Functions/Test_SpecialFunctions.cs b/Test_Numerics/Mathematics/Special Functions/Test_SpecialFunctions.cs index b5372e34..9ca6b01d 100644 --- a/Test_Numerics/Mathematics/Special Functions/Test_SpecialFunctions.cs +++ b/Test_Numerics/Mathematics/Special Functions/Test_SpecialFunctions.cs @@ -1,4 +1,4 @@ -using System; +using System; using Microsoft.VisualStudio.TestTools.UnitTesting; using Numerics.Mathematics.SpecialFunctions; @@ -222,7 +222,7 @@ public void Test_CombinationsNum() public void Test_Combinations() { int[,] valid = new int[7,3] - { + { { 1, 0 ,0 }, { 0, 1, 0 }, { 0, 0, 1 }, @@ -307,4 +307,4 @@ public void Test_PolynomialRev_1() Assert.AreEqual(valid, actual); } } -} \ No newline at end of file +} diff --git a/Test_Numerics/Sampling/MCMC/Test_Gibbs.cs b/Test_Numerics/Sampling/MCMC/Test_Gibbs.cs index 028a2785..a100798f 100644 --- a/Test_Numerics/Sampling/MCMC/Test_Gibbs.cs +++ b/Test_Numerics/Sampling/MCMC/Test_Gibbs.cs @@ -36,10 +36,16 @@ public void Test_Gibbs_NormalDist_RStan() int n = sample.Length; var mu = Statistics.Mean(sample); double mu0 = 0, sigma0 = 5E5; - var muPrior = new Normal(mu0, sigma0); - double alpha0 = 2, beta0 = 0.001; - var sigmaPrior = new InverseGamma(beta0, alpha0); - var priors = new List { muPrior, sigmaPrior }; + // Preserve the reference model's limiting non-informative inverse-gamma prior. + double variancePriorShape = 0d, variancePriorScale = 0d; + + // Use proper distributions solely to initialize the sampler state. + double initializationShape = 2d, initializationScale = 0.001d; + var muInitializationPrior = new Normal(mu0, sigma0); + var sigmaInitializationPrior = new InverseGamma(initializationScale, initializationShape); + var conditionalMean = new Normal(); + var conditionalVariance = new InverseGamma(); + var priors = new List { muInitializationPrior, sigmaInitializationPrior }; // Create log-likelihood function double logLH(double[] x) @@ -51,20 +57,15 @@ double logLH(double[] x) // Create proposal function double[] proposal(double[] x, Random random) { - // Sample mu - double mun = (n * mu + mu0 / 2) / (n + 1 / (sigma0 * sigma0)); - double sigma2 = (x[1] * x[1]) / (n + (x[1] * x[1]) / (sigma0 * sigma0)); - muPrior.SetParameters(mun, Math.Sqrt(sigma2)); - double mup = muPrior.InverseCDF(random.NextDouble()); + // Sample the conditional mean given the current standard deviation. + var meanParameters = ConditionalMeanParameters(n, mu, x[1], mu0, sigma0); + conditionalMean.SetParameters(meanParameters.Mean, meanParameters.StandardDeviation); + double mup = conditionalMean.InverseCDF(random.NextDouble()); - // Sample sigma - double alpha1 = n / 2d; - double sse = 0; - for (int i = 0; i < sample.Length; i++) - sse += Math.Pow(sample[i] - mup, 2); - double beta1 = sse / 2d; - sigmaPrior.SetParameters(new double[] { beta1, alpha1 }); - double sig2p = sigmaPrior.InverseCDF(random.NextDouble()); + // Sample the conditional variance, then return its square root as sigma. + var varianceParameters = ConditionalVarianceParameters(sample, mup, variancePriorShape, variancePriorScale); + conditionalVariance.SetParameters(new[] { varianceParameters.Scale, varianceParameters.Shape }); + double sig2p = conditionalVariance.InverseCDF(random.NextDouble()); // return proposal vector return new double[] { mup, Math.Sqrt(sig2p) }; @@ -74,6 +75,10 @@ double[] proposal(double[] x, Random random) var sampler = new Gibbs(priors, logLH, proposal); sampler.Sample(); var results = new MCMCResults(sampler); + Assert.AreEqual(mu0, muInitializationPrior.Mu); + Assert.AreEqual(sigma0, muInitializationPrior.Sigma); + Assert.AreEqual(initializationScale, sigmaInitializationPrior.Beta); + Assert.AreEqual(initializationShape, sigmaInitializationPrior.Alpha); /* Below are the results from 'rstan' using comparable MCMC settings: * mean se_mean sd 5% 50% 95% n_eff Rhat @@ -82,8 +87,8 @@ double[] proposal(double[] x, Random random) * lp__ -466.13 0.01 1.03 -468.17 -465.81 -465.15 9958 1 * * Since MCMC methods rely on random number generation, results will not be - * exactly the same as those produced by other samplers. Therefore, these - * comparisons aim to verify whether the results are within 5% of 'rstan' results. + * exactly the same as those produced by other samplers. Therefore, these + * comparisons verify that the results remain within 5% of the 'rstan' results. */ // Mu @@ -99,6 +104,69 @@ double[] proposal(double[] x, Random random) Assert.AreEqual(4796.63, results.ParameterResults[1].SummaryStatistics.Median, 0.05 * 4796.63); Assert.AreEqual(5771.81, results.ParameterResults[1].SummaryStatistics.UpperCI, 0.05 * 5771.81); } + + /// + /// Verifies the Normal and inverse-gamma full-conditionals with an informative, nonzero prior. + /// + [TestMethod] + public void Test_ConditionalParameters_InformativePrior() + { + var meanParameters = ConditionalMeanParameters(4, 10d, 2d, 2d, 3d); + Assert.AreEqual(9.2d, meanParameters.Mean, 1E-12); + Assert.AreEqual(Math.Sqrt(0.9d), meanParameters.StandardDeviation, 1E-12); + + var varianceParameters = ConditionalVarianceParameters(new[] { 8d, 10d, 12d, 14d }, 10d, 2d, 0.5d); + Assert.AreEqual(12.5d, varianceParameters.Scale, 1E-12); + Assert.AreEqual(4d, varianceParameters.Shape, 1E-12); + } + + /// + /// Computes the Normal full-conditional parameters for the population mean when the + /// likelihood variance is fixed at its current Gibbs state. + /// + /// The number of observations. + /// The arithmetic mean of the observations. + /// The current likelihood standard deviation. + /// The Normal prior mean. + /// The Normal prior standard deviation. + /// The conditional Normal mean and standard deviation. + private static (double Mean, double StandardDeviation) ConditionalMeanParameters( + int sampleSize, + double sampleMean, + double currentStandardDeviation, + double priorMean, + double priorStandardDeviation) + { + double likelihoodVariance = currentStandardDeviation * currentStandardDeviation; + double priorVariance = priorStandardDeviation * priorStandardDeviation; + double posteriorVariance = 1d / (sampleSize / likelihoodVariance + 1d / priorVariance); + double posteriorMean = posteriorVariance * + (sampleSize * sampleMean / likelihoodVariance + priorMean / priorVariance); + return (posteriorMean, Math.Sqrt(posteriorVariance)); + } + + /// + /// Computes the inverse-gamma full-conditional parameters for the likelihood variance. + /// + /// The observed sample. + /// The mean drawn during the current Gibbs iteration. + /// The inverse-gamma prior shape. + /// The inverse-gamma prior scale. + /// The conditional inverse-gamma scale and shape, in the order expected by . + private static (double Scale, double Shape) ConditionalVarianceParameters( + IList sample, + double conditionalMean, + double priorShape, + double priorScale) + { + double sumOfSquaredErrors = 0d; + for (int i = 0; i < sample.Count; i++) + { + double residual = sample[i] - conditionalMean; + sumOfSquaredErrors += residual * residual; + } + return (priorScale + sumOfSquaredErrors / 2d, priorShape + sample.Count / 2d); + } } } diff --git a/Test_Numerics/Sampling/MCMC/Test_MCMCInitialization.cs b/Test_Numerics/Sampling/MCMC/Test_MCMCInitialization.cs new file mode 100644 index 00000000..802d8204 --- /dev/null +++ b/Test_Numerics/Sampling/MCMC/Test_MCMCInitialization.cs @@ -0,0 +1,178 @@ +using System; +using System.Collections.Generic; +using Microsoft.VisualStudio.TestTools.UnitTesting; +using Numerics.Distributions; +using Numerics.Mathematics.LinearAlgebra; +using Numerics.Mathematics.Optimization; +using Numerics.Sampling.MCMC; + +namespace Sampling.MCMC +{ + /// + /// Regression tests for MCMC initialization behavior. + /// + [TestClass] + public class Test_MCMCInitialization + { + /// + /// Verifies that MAP initialization stores the original log-likelihood sign without an extra evaluation. + /// + [TestMethod] + public void MapInitializationStoresLogLikelihoodWithoutReevaluation() + { + int optimizerEvaluations = 0; + var optimizer = new DifferentialEvolution( + parameters => + { + optimizerEvaluations++; + return QuadraticLogLikelihood(parameters); + }, + 1, + new[] { -5d }, + new[] { 5d }) + { + ReportFailure = false + }; + optimizer.Maximize(); + Assert.AreEqual(OptimizationStatus.Success, optimizer.Status); + + int samplerEvaluations = 0; + var sampler = new InitializableRwmh( + new List { new Uniform(-5d, 5d) }, + parameters => + { + samplerEvaluations++; + return QuadraticLogLikelihood(parameters); + }) + { + NumberOfChains = 1, + InitialIterations = 2, + Initialize = MCMCSampler.InitializationType.MAP + }; + + sampler.InitializeOnly(); + + Assert.IsFalse(sampler.MAPInitializationFailed); + Assert.AreEqual(optimizerEvaluations + sampler.InitialIterations, samplerEvaluations); + Assert.AreEqual(QuadraticLogLikelihood(sampler.MAP.Values), sampler.MAP.Fitness, 1E-12); + Assert.AreEqual(-3d, sampler.MAP.Fitness, 1E-6); + } + + /// + /// Verifies that NUTS uses an analytical gradient during its step-size initialization heuristic. + /// + [TestMethod] + public void NutsInitializationUsesConfiguredGradientAndReducesLikelihoodWork() + { + var priors = new List { new Uniform(-5d, 5d) }; + int analyticalLikelihoodEvaluations = 0; + int analyticalGradientEvaluations = 0; + var analytical = new InitializableNuts( + priors, + parameters => + { + analyticalLikelihoodEvaluations++; + return -parameters[0] * parameters[0]; + }, + parameters => + { + analyticalGradientEvaluations++; + return new Vector(new[] { -2d * parameters[0] }); + }); + ConfigureSingleChainInitialization(analytical); + analytical.InitializeOnly(); + + int numericalLikelihoodEvaluations = 0; + var numerical = new InitializableNuts( + priors, + parameters => + { + numericalLikelihoodEvaluations++; + return -parameters[0] * parameters[0]; + }); + ConfigureSingleChainInitialization(numerical); + numerical.InitializeOnly(); + + Assert.IsGreaterThan(0, analyticalGradientEvaluations); + Assert.IsGreaterThan(analyticalLikelihoodEvaluations, numericalLikelihoodEvaluations); + } + + /// + /// Evaluates the quadratic log-likelihood used to verify MAP fitness and evaluation counts. + /// + /// The parameter vector to evaluate. + /// The quadratic log-likelihood centered at one with a maximum of negative three. + private static double QuadraticLogLikelihood(double[] parameters) + { + double difference = parameters[0] - 1d; + return -3d - difference * difference; + } + + /// + /// Applies deterministic, minimal settings to an initialization-only sampler. + /// + /// The sampler to configure. + private static void ConfigureSingleChainInitialization(MCMCSampler sampler) + { + sampler.NumberOfChains = 1; + sampler.InitialIterations = 1; + sampler.PRNGSeed = 12345; + } + + /// + /// Exposes random-walk Metropolis-Hastings initialization for focused regression testing. + /// + private sealed class InitializableRwmh : RWMH + { + /// + /// Initializes a new instance of the class. + /// + /// The prior distributions for the sampled parameters. + /// The log-likelihood function to evaluate. + internal InitializableRwmh( + List priorDistributions, + LogLikelihood logLikelihood) + : base(priorDistributions, logLikelihood, new Matrix(1)) + { + } + + /// + /// Initializes the chains without running the sampler. + /// + /// The initialized chain parameter sets. + internal ParameterSet[] InitializeOnly() + { + return InitializeChains(); + } + } + + /// + /// Exposes No-U-Turn sampler initialization for focused gradient-routing tests. + /// + private sealed class InitializableNuts : NUTS + { + /// + /// Initializes a new instance of the class. + /// + /// The prior distributions for the sampled parameters. + /// The log-likelihood function to evaluate. + /// The optional analytical gradient of the log-likelihood. + internal InitializableNuts( + List priorDistributions, + LogLikelihood logLikelihood, + HMC.Gradient gradient = null) + : base(priorDistributions, logLikelihood, maxTreeDepth: 2, gradientFunction: gradient) + { + } + + /// + /// Initializes the chains and sampler-specific step-size state without drawing samples. + /// + internal void InitializeOnly() + { + _chainStates = InitializeChains(); + InitializeCustomSettings(); + } + } + } +} diff --git a/Test_Numerics/Sampling/MCMC/Test_MCMCResults_Recompute.cs b/Test_Numerics/Sampling/MCMC/Test_MCMCResults_Recompute.cs index 5a8ec980..b0c0fc54 100644 --- a/Test_Numerics/Sampling/MCMC/Test_MCMCResults_Recompute.cs +++ b/Test_Numerics/Sampling/MCMC/Test_MCMCResults_Recompute.cs @@ -80,6 +80,9 @@ private static (double[] rhats, double[] esss, double[][,] acfs) SeedDiagnostics return (rhats, esss, acfs); } + /// + /// Verifies recomputation preserves the output sample reference. + /// [TestMethod] public void Recompute_PreservesOutputReference() { @@ -93,6 +96,9 @@ public void Recompute_PreservesOutputReference() Assert.HasCount(1000, results.Output, "Output count must be preserved."); } + /// + /// Verifies recomputation preserves MAP results. + /// [TestMethod] public void Recompute_PreservesMAP() { @@ -108,6 +114,9 @@ public void Recompute_PreservesMAP() Assert.AreEqual(fitnessBefore, results.MAP.Fitness, 1e-12, "MAP fitness must not change."); } + /// + /// Verifies recomputation preserves R-hat diagnostics. + /// [TestMethod] public void Recompute_PreservesRhat() { @@ -123,6 +132,9 @@ public void Recompute_PreservesRhat() } } + /// + /// Verifies recomputation preserves effective sample size diagnostics. + /// [TestMethod] public void Recompute_PreservesESS() { @@ -138,6 +150,9 @@ public void Recompute_PreservesESS() } } + /// + /// Verifies recomputation preserves autocorrelation diagnostics. + /// [TestMethod] public void Recompute_PreservesAutocorrelation() { @@ -162,6 +177,9 @@ public void Recompute_PreservesAutocorrelation() } } + /// + /// Verifies recomputation narrows credible intervals when alpha increases. + /// [TestMethod] public void Recompute_NarrowsCIWhenAlphaIncreases() { @@ -183,6 +201,9 @@ public void Recompute_NarrowsCIWhenAlphaIncreases() } } + /// + /// Verifies recomputation widens credible intervals when alpha decreases. + /// [TestMethod] public void Recompute_WidensCIWhenAlphaDecreases() { @@ -204,6 +225,9 @@ public void Recompute_WidensCIWhenAlphaDecreases() } } + /// + /// Verifies recomputation preserves mean and median summaries. + /// [TestMethod] public void Recompute_PreservesMeanAndMedian() { @@ -223,6 +247,9 @@ public void Recompute_PreservesMeanAndMedian() } } + /// + /// Verifies recomputation on empty results does not throw. + /// [TestMethod] public void Recompute_OnEmptyResults_DoesNotThrow() { diff --git a/Test_Numerics/Utilities/Test_Tools.cs b/Test_Numerics/Utilities/Test_Tools.cs index 1287b891..59dc4a33 100644 --- a/Test_Numerics/Utilities/Test_Tools.cs +++ b/Test_Numerics/Utilities/Test_Tools.cs @@ -1,6 +1,9 @@ using Microsoft.VisualStudio.TestTools.UnitTesting; using Numerics; +using System; using System.Collections.Generic; +using System.Threading; +using System.Threading.Tasks; namespace Utilities { @@ -335,6 +338,155 @@ public void Test_MaxIndicator() Assert.AreEqual(1, result); } + /// + /// Verifies concurrent additions are committed exactly once and return their serialization positions. + /// + [TestMethod] + public void Test_ParallelAdd_CommitsEveryConcurrentUpdate() + { + const int updateCount = 25000; + double total = 0d; + var committedValues = new double[updateCount]; + + Parallel.For( + 0, + updateCount, + index => committedValues[index] = Tools.ParallelAdd(ref total, 1d)); + + Assert.AreEqual((double)updateCount, total); + Array.Sort(committedValues); + for (int i = 0; i < committedValues.Length; i++) + { + Assert.AreEqual(i + 1d, committedValues[i]); + } + } + + /// + /// Verifies the compare-and-swap loop terminates when the accumulator contains NaN. + /// + [TestMethod] + public void Test_ParallelAdd_NaNAccumulatorTerminates() + { + double total = double.NaN; + + double result = Tools.ParallelAdd(ref total, 1d); + + Assert.IsTrue(double.IsNaN(total)); + Assert.IsTrue(double.IsNaN(result)); + } + + /// + /// Verifies signed zero, NaN, and infinity follow double-precision addition semantics. + /// + [TestMethod] + public void Test_ParallelAdd_HandlesSpecialValues() + { + double negativeZero = BitConverter.Int64BitsToDouble(long.MinValue); + + double total = negativeZero; + double result = Tools.ParallelAdd(ref total, negativeZero); + + Assert.AreEqual(long.MinValue, BitConverter.DoubleToInt64Bits(total)); + Assert.AreEqual(long.MinValue, BitConverter.DoubleToInt64Bits(result)); + + total = negativeZero; + result = Tools.ParallelAdd(ref total, 0d); + + Assert.AreEqual(0L, BitConverter.DoubleToInt64Bits(total)); + Assert.AreEqual(0L, BitConverter.DoubleToInt64Bits(result)); + + total = 1d; + result = Tools.ParallelAdd(ref total, double.NaN); + + Assert.IsTrue(double.IsNaN(total)); + Assert.IsTrue(double.IsNaN(result)); + + total = double.PositiveInfinity; + result = Tools.ParallelAdd(ref total, 1d); + + Assert.AreEqual(double.PositiveInfinity, total); + Assert.AreEqual(double.PositiveInfinity, result); + + total = double.NegativeInfinity; + result = Tools.ParallelAdd(ref total, -1d); + + Assert.AreEqual(double.NegativeInfinity, total); + Assert.AreEqual(double.NegativeInfinity, result); + + total = double.PositiveInfinity; + result = Tools.ParallelAdd(ref total, double.NegativeInfinity); + + Assert.IsTrue(double.IsNaN(total)); + Assert.IsTrue(double.IsNaN(result)); + } + + /// + /// Verifies concurrent mixed-sign additions are all committed when the accumulator revisits values. + /// + [TestMethod] + public void Test_ParallelAdd_CommitsMixedSignUpdates() + { + const int updateCount = 25000; + double total = 0d; + + Parallel.For( + 0, + updateCount, + index => Tools.ParallelAdd(ref total, index % 2 == 0 ? 2d : -1d)); + + Assert.AreEqual(updateCount / 2d, total); + } + + /// + /// Verifies a negative-to-positive zero race cannot cause finite concurrent updates to be lost. + /// + [TestMethod] + public void Test_ParallelAdd_SignedZeroContentionDoesNotLoseUpdates() + { + const int rounds = 256; + double total = 0d; + bool everyUpdateCommitted = true; + double negativeZero = BitConverter.Int64BitsToDouble(long.MinValue); + + using (var startBarrier = new Barrier(5)) + using (var finishBarrier = new Barrier(5)) + { + Func createWorker = valueToAdd => Task.Factory.StartNew( + () => + { + for (int round = 0; round < rounds; round++) + { + startBarrier.SignalAndWait(); + Tools.ParallelAdd(ref total, valueToAdd); + finishBarrier.SignalAndWait(); + } + }, + CancellationToken.None, + TaskCreationOptions.LongRunning, + TaskScheduler.Default); + + Task[] workers = + { + createWorker(0d), + createWorker(0d), + createWorker(1d), + createWorker(1d) + }; + + for (int round = 0; round < rounds; round++) + { + total = negativeZero; + startBarrier.SignalAndWait(); + finishBarrier.SignalAndWait(); + if (total != 2d) everyUpdateCommitted = false; + } + + Task.WaitAll(workers); + } + + Assert.IsTrue(everyUpdateCommitted, "A finite update was lost during signed-zero contention."); + } + /// /// Testing the Log-sum-exponential function. /// diff --git a/codemeta.json b/codemeta.json index 046bf4b6..a3be965b 100644 --- a/codemeta.json +++ b/codemeta.json @@ -4,9 +4,9 @@ "name": "Numerics: A .NET Library for Numerical Computing, Statistical Analysis, and Risk Assessment", "alternateName": "Numerics", "description": "A free and open-source .NET library providing numerical methods, probability distributions, statistical analysis, and Bayesian inference tools for quantitative risk assessment in water resources engineering.", - "version": "2.1.2", + "version": "2.1.3", "dateCreated": "2023-09-28", - "dateModified": "2026-07-09", + "dateModified": "2026-07-15", "license": "https://spdx.org/licenses/0BSD", "codeRepository": "https://github.com/USACE-RMC/Numerics", "issueTracker": "https://github.com/USACE-RMC/Numerics/issues",