Fix test regressions and refactor NoncentralT.ZBrent to use Brent.Solve

- Revert RSR denominator from N-1 back to N (ensures RSR=1.0 for
  constant-at-mean predictions); fix misleading comment
- Fix assertion argument order in Test_InverseChiSquared,
  Test_TimeSeries, Test_kNN, and Test_MCMCDiagnostics
- Replace 120-line private ZBrent method in NoncentralT with a single
  call to Brent.Solve, using Math.Min/Max for bracket ordering

Author: HadenSmith

Numerics/Data/Statistics/GoodnessOfFit.cs

@@ -804,8 +804,8 @@ public static double PBIAS(IList<double> observedValues, IList<double> modeledVa
         /// Lower RSR values indicate better model performance.
         /// </para>
         /// <para>
-        /// <b>Note:</b> This method uses sample standard deviation (N-1 denominator) for the observed values,
-        /// which is consistent with the R hydroGOF package implementation.
+        /// <b>Note:</b> This method uses the population standard deviation (N denominator) for both RMSE and
+        /// observed variability, ensuring RSR = 1.0 when predictions equal the observed mean.
         /// </para>
         /// </remarks>
         public static double RSR(IList<double> observedValues, IList<double> modeledValues)
@@ -828,13 +828,13 @@ public static double RSR(IList<double> observedValues, IList<double> modeledValu
             double rmse = Math.Sqrt(sse / n);
             double obsMean = sumObs / n;
 
-            // Second pass to compute sample standard deviation (using N-1)
+            // Second pass to compute standard deviation
             double variance = 0d;
             for (int i = 0; i < n; i++)
             {
                 variance += Tools.Sqr(observedValues[i] - obsMean);
             }
-            double stdDev = Math.Sqrt(variance / (n - 1));
+            double stdDev = Math.Sqrt(variance / n);
 
             return rmse / stdDev;
         }

Numerics/Distributions/Univariate/NoncentralT.cs

@@ -31,6 +31,7 @@
 using System;
 using System.Collections.Generic;
 using Numerics.Mathematics.Optimization;
+using Numerics.Mathematics.RootFinding;
 using Numerics.Mathematics.SpecialFunctions;
 
 namespace Numerics.Distributions
@@ -566,7 +567,7 @@ private double NCT_INV(double p, double df, double delta)
                 iter = iter + 1;
             }
             // Solve for T using Brent
-            double ANS = ZBrent(t0, t1, y0, y1, xtol, df, delta, p);
+            double ANS = Brent.Solve(x => NCTDist(x, df, delta) - p, Math.Min(t0, t1), Math.Max(t0, t1), xtol, reportFailure: false);
             return ANS;
         }
 
@@ -597,128 +598,6 @@ private double NCTInv0(double P, double N, double D)
             }
         }
 
-        private double ZBrent(double X1, double X2, double y1, double y2, double Tol, double N, double Dnc, double Perc)
-        {
-            double ZBrentRet = default;
-            // 
-            // Finds the zero of NCTDist(x, n, d) - p given that [x1, x2] brackets the zero and
-            // Y1 = value at X1, Y2 = value at X2.
-            // 
-            // Translated from Numerical Recipes (1986).
-            // William A. Huber, 24 March 2001.
-            // 
-            double a;
-            double b;
-            var c = default(double);
-            double fc;
-            var D = default(double);
-            var e = default(double);
-            double tol1;
-            double xm;
-            double s;
-            double P;
-            double q;
-            double r;
-            double fa;
-            double fb;
-            int iter;
-            const int itmax = 100;
-            const double eps = 0.00000003d;
-            a = X1;
-            b = X2;
-            fa = y1;
-            fb = y2;
-            if (fb * fa > 0d)
-            {
-                throw new ArgumentException("Brent's method failed because the root is not bracketed.");
-            }
-
-            fc = fb;
-            for (iter = 1; iter <= itmax; iter++)
-            {
-                if (fb * fc > 0d)
-                {
-                    c = a;
-                    fc = fa;
-                    D = b - a;
-                    e = D;
-                }
-
-                if (Math.Abs(fc) < Math.Abs(fb))
-                {
-                    a = b;
-                    b = c;
-                    c = a;
-                    fa = fb;
-                    fb = fc;
-                    fc = fa;
-                }
-
-                tol1 = 2.0d * eps * Math.Abs(b) + 0.5d * Tol;
-                xm = 0.5d * (c - b);
-                if (Math.Abs(xm) <= tol1 || fb == 0d)
-                {
-                    return b;
-                }
-
-                if (Math.Abs(e) >= tol1 && Math.Abs(fa) > Math.Abs(fb))
-                {
-                    s = fb / fa;
-                    if (a == c)
-                    {
-                        P = 2.0d * xm * s;
-                        q = 1.0d - s;
-                    }
-                    else
-                    {
-                        q = fa / fc;
-                        r = fb / fc;
-                        P = s * (2.0d * xm * q * (q - r) - (b - a) * (r - 1.0d));
-                        q = (q - 1.0d) * (r - 1.0d) * (s - 1.0d);
-                    }
-
-                    if (P > 0d)
-                        q = -q;
-                    P = Math.Abs(P);
-                    if (2.0d * P < 3.0d * xm * q - Math.Abs(tol1 * q) & 2.0d * P < Math.Abs(e * q))
-                    {
-                        e = D;
-                        D = P / q;
-                    }
-                    else
-                    {
-                        D = xm;
-                        e = D;
-                    }
-                }
-                else
-                {
-                    D = xm;
-                    e = D;
-                }
-
-                a = b;
-                fa = fb;
-                if (Math.Abs(D) > tol1)
-                {
-                    b = b + D;
-                }
-                else if (xm > 0d)  // b = b + SIGN(tol1, xm)
-                {
-                    b = b + Math.Abs(tol1);
-                }
-                else
-                {
-                    b = b - Math.Abs(tol1);
-                }
-
-                fb = NCTDist(b, N, Dnc) - Perc;
-            }
-
-            ZBrentRet = b;
-            return ZBrentRet;
-        }
-
         /// <inheritdoc/>
         public override UnivariateDistributionBase Clone()
         {

Test_Numerics/Data/Time Series/Test_TimeSeries.cs

@@ -1008,8 +1008,8 @@ public void Test_SeasonalDecompose()
             var (trend, seasonal, residual) = ts.SeasonalDecompose(period);
 
             // Verify trend is returned
-            Assert.IsGreaterThan(trend.Count, 0, "Trend should have values.");
-            Assert.IsLessThanOrEqualTo(trend.Count, n, "Trend should not exceed original length.");
+            Assert.IsGreaterThan(0, trend.Count, "Trend should have values.");
+            Assert.IsLessThanOrEqualTo(n, trend.Count, "Trend should not exceed original length.");
 
             // Verify seasonal has correct length
             // HasCount does not accept double[] (not IEnumerable), so suppress MSTEST0037
@@ -1018,7 +1018,7 @@ public void Test_SeasonalDecompose()
 #pragma warning restore MSTEST0037
 
             // Verify residual is returned
-            Assert.IsGreaterThan(residual.Count, 0, "Residual should have values.");
+            Assert.IsGreaterThan(0, residual.Count, "Residual should have values.");
 
             // Verify the decomposition is additive: original = trend + seasonal + residual
             // This is an exact mathematical identity by construction
@@ -1043,7 +1043,7 @@ public void Test_SeasonalDecompose()
             // Verify seasonal values are not all zero (they should have nonzero amplitude)
             double maxSeasonal = seasonal.Max();
             double minSeasonal = seasonal.Min();
-            Assert.IsGreaterThan(maxSeasonal - minSeasonal, 1.0,
+            Assert.IsGreaterThan(1.0, maxSeasonal - minSeasonal,
                 "Seasonal component should have nonzero amplitude for sinusoidal input.");
         }
 

Test_Numerics/Distributions/Univariate/Test_InverseChiSquared.cs

@@ -266,8 +266,8 @@ public void Test_CDF_InverseCDF_Roundtrip()
             Assert.AreEqual(0.8911780189, dist.CDF(2.0), 1E-6);
 
             // CDF must be non-decreasing
-            Assert.IsLessThan(dist.CDF(0.5), dist.CDF(1.0));
-            Assert.IsLessThan(dist.CDF(1.0), dist.CDF(2.0));
+            Assert.IsLessThan(dist.CDF(1.0), dist.CDF(0.5));
+            Assert.IsLessThan(dist.CDF(2.0), dist.CDF(1.0));
 
             // InverseCDF values from v*sigma / scipy.stats.chi2.isf(p, v)
             Assert.AreEqual(0.6255012152, dist.InverseCDF(0.1), 1E-4);

Test_Numerics/Machine Learning/Supervised/Test_kNN.cs

@@ -190,12 +190,12 @@ public void Test_GetNeighbors_MultiRow()
             Assert.IsNotNull(neighbors);
 
             // First query (0.5, 0.5): nearest neighbors should be from cluster A (indices 0-5)
-            Assert.IsLessThan(neighbors[0], 6, $"First query's nearest neighbor should be in cluster A, got index {neighbors[0]}");
-            Assert.IsLessThan(neighbors[1], 6, $"First query's 2nd nearest should be in cluster A, got index {neighbors[1]}");
+            Assert.IsLessThan(6, neighbors[0], $"First query's nearest neighbor should be in cluster A, got index {neighbors[0]}");
+            Assert.IsLessThan(6, neighbors[1], $"First query's 2nd nearest should be in cluster A, got index {neighbors[1]}");
 
             // Second query (100.5, 100.5): nearest neighbors should be from cluster B (indices 6-11)
-            Assert.IsGreaterThanOrEqualTo(neighbors[2], 6, $"Second query's nearest neighbor should be in cluster B, got index {neighbors[2]}");
-            Assert.IsGreaterThanOrEqualTo(neighbors[3], 6, $"Second query's 2nd nearest should be in cluster B, got index {neighbors[3]}");
+            Assert.IsGreaterThanOrEqualTo(6, neighbors[2], $"Second query's nearest neighbor should be in cluster B, got index {neighbors[2]}");
+            Assert.IsGreaterThanOrEqualTo(6, neighbors[3], $"Second query's 2nd nearest should be in cluster B, got index {neighbors[3]}");
         }
 
     }

Test_Numerics/Sampling/MCMC/Test_MCMCDiagnostics.cs

@@ -76,14 +76,14 @@ public void Test_GelmanRubin_WithWarmup()
 
             // Without warmup, R-hat should be high (chains have different initial distributions)
             var rhatNoWarmup = MCMCDiagnostics.GelmanRubin(chains, 0);
-            Assert.IsGreaterThan(rhatNoWarmup[0], 1.1, $"R-hat without warmup should be > 1.1, got {rhatNoWarmup[0]}");
+            Assert.IsGreaterThan(1.1, rhatNoWarmup[0], $"R-hat without warmup should be > 1.1, got {rhatNoWarmup[0]}");
 
             // With warmup=50, R-hat should be close to 1.0 (converged portion only)
             var rhatWithWarmup = MCMCDiagnostics.GelmanRubin(chains, 50);
-            Assert.IsLessThan(rhatWithWarmup[0], 1.1, $"R-hat with warmup=50 should be < 1.1, got {rhatWithWarmup[0]}");
+            Assert.IsLessThan(1.1, rhatWithWarmup[0], $"R-hat with warmup=50 should be < 1.1, got {rhatWithWarmup[0]}");
 
             // Warmup should improve R-hat (make it closer to 1.0)
-            Assert.IsLessThan(rhatWithWarmup[0], rhatNoWarmup[0],
+            Assert.IsLessThan(rhatNoWarmup[0], rhatWithWarmup[0],
                 "R-hat with warmup should be closer to 1.0 than without warmup");
         }