Comprehensive bugfixes, correctness improvements, and code quality enhancements

Phase 1 - Critical correctness fixes:
- Fix InverseChiSquared CDF and InverseCDF formulas (was computing survival function)
- Fix NormalCopula CDF computing survival copula instead of copula
- Fix AMH, Frank, and Normal copula ValidateParameter returning non-null for valid params
- Fix MultipleGrubbsBeckTest static mutable fields (thread safety)
- Fix TimeSeries OneQuarter interval off by factor of 4
- Fix Search.Hunt always assuming descending order for OrderedPairedData overload
- Fix Bilinear UseSmartSearch copy-paste bug (X2LI never set)
- Fix KMeans++ initialization always using centroid[0]
- Fix SNIS weight filter wrong boolean operator (|| vs &&)
- Fix KNearestNeighbors multi-row GetNeighbors overwriting results

Phase 2 - High-priority fixes:
- Fix TimeSeries StartDate/EndDate crash on empty series
- Fix TimeSeries InterpolateMissingData index -1 crash
- Fix TimeSeries CumulativeSum NaN vs zero confusion
- Fix TimeSeries Divide silent infinity on zero divisor
- Fix TimeSeries StandardDeviation/Duration/SummaryPercentiles NaN handling
- Fix LogPearsonTypeIII Clone() dropping Base parameter
- Fix ChiSquared PDF overflow for large degrees of freedom (log-space)
- Fix Bernoulli Skewness/Kurtosis division by zero at p=0,1
- Fix Binomial InverseCDF off-by-one (never returns NumberOfTrials)
- Fix Rayleigh MoM wrong estimator
- Fix UniformDiscrete StandardDeviation wrong formula (continuous vs discrete)
- Fix GeneralizedBeta Mode 0/0 at Alpha=Beta=1
- Fix TruncatedNormal ValidateParameters reading property instead of parameter
- Fix EmpiricalDistribution/KernelDensity GetParameters NotImplementedException
- Fix Statistics HarmonicMean missing zero guard
- Fix GoodnessOfFit RSR using N instead of N-1
- Fix Histogram multiple issues (BinWidth=0, swapped args, empty guard)
- Fix MCMCDiagnostics GelmanRubin wrong denominator with warmup
- Fix GLM bounds inversion for negative intercept
- Fix HMC boundary enforcement (add momentum reflection)
- Fix NaiveBayes NaN for singleton classes
- Fix JenksCluster numerically unstable variance (Welford's algorithm)
- Fix ARWMH ProposalSigma division by zero when N=1
- Fix PowerFunction MinimumOfParameters/MaximumOfParameters wrong array length
- Fix LinearFunction/PowerFunction SetParameters not updating _normal
- Fix LinearFunction InverseFunction division by zero when Beta=0
- Fix MCMCSampler AcceptanceRates/ProgressChanged division by zero
- Fix DEMCz infinite loop when PopulationMatrix.Count=1
- Fix Autocorrelation/Correlation division by zero for constant series
- Fix HypothesisTests division-by-zero issues
- Fix BoxCox/GrubbsBeckTest non-positive values not guarded
- Remove VonMises Kurtosis dead computation

Phase 3 - Security (TimeSeriesDownload):
- Use static HttpClient instances (prevent socket exhaustion)
- Add CancellationToken support to all public async methods
- Validate redirect URLs against expected domains
- Remove UseProxy=false bypass
- Replace browser User-Agent with library identifier
- Replace Google.com connectivity check with target endpoint
- Sanitize GHCN station numbers
- Fix temp file leak on error
- Add explicit timeouts and retry logic with exponential backoff
- Use automatic gzip decompression

Phase 4 - Robustness:
- Fix Brent.Solve returning 0 instead of last iterate on failure
- Fix Series.SyncRoot returning unstable object
- Fix FrankCopula generator division by zero at Theta=0
- Fix TruncatedDistribution ValidateParameters exact float equality
- Fix Mixture EM convergence check placement
- Fix Deterministic/Binomial wrong parameter names in exceptions
- Fix JsonConverters silent zero-fill on malformed input
- Fix SafeProgressReporter MessageCount thread safety
- Fix ExtensionMethods.NextDoubles exact type check (use 'is')
- Fix TimeSeries FillMissingDates O(n^2) to O(n)
- Fix TimeSeries ResampleWithKNN index exceeding neighbor count
- Fix AdaptiveGaussKronrod wrong exception type

Phase 5 - Code quality:
- Remove useless catch/throw blocks (SobolSequence)
- Remove unused Microsoft.VisualBasic using (Gamma.cs)
- Fix CubicSpline XML doc saying "linear interpolation"
- Remove dead code (HypothesisTests, Search, YeoJohnson, DEMCzs)
- Fix KernelDensity typo GuassianKernel -> GaussianKernel
- Fix bitwise | vs logical || (Ordinate, Search)
- Remove unused NUnit references from test project
- Fix KNN XML doc "clusters" -> "nearest neighbors"
- Fix PowerFunction.Minimum setter no-op -> NotSupportedException
- Enable CS1591 (missing XML doc) warning and add missing XML docs
- Add XML docs for BinaryHeap, Network, SafeProgressReporter,
  MCMCSampler, MultivariateNormal, KernelDensity, SimulatedAnnealing,
  TimeSeries.FillMissingDates

Phase 6 - Unit tests with Python-verified reference values:
- Update InverseChiSquared tests for corrected CDF/InverseCDF
- Add InverseChiSquared CDF/InverseCDF roundtrip test
- Add NormalCopula CDF asymmetric inputs test
- Add Bernoulli boundary skewness/kurtosis test
- Add Binomial InverseCDF boundary test
- Add ChiSquared PDF large DoF test
- Add GeneralizedBeta Mode uniform case test
- Add Rayleigh MoM estimation test
- Add UniformDiscrete StandardDeviation test
- Add Statistics HarmonicMean zero value test
- Add KNN multi-row GetNeighbors test
- Add MCMCDiagnostics GelmanRubin warmup test
- Add HMC non-finite gradient crash test
- Fix MSTEST0037 analyzer warnings across test project

Author: HadenSmith

Numerics/Data/Interpolation/Bilinear.cs

@@ -119,7 +119,7 @@ public bool UseSmartSearch
             {
                 _useSmartSearch = value;
                 X1LI.UseSmartSearch = value;
-                X1LI.UseSmartSearch = value;
+                X2LI.UseSmartSearch = value;
             }
         }
 

Numerics/Data/Interpolation/CubicSpline.cs

@@ -62,7 +62,7 @@ public class CubicSpline : Interpolater
     {
 
         /// <summary>
-        /// Construct new linear interpolation. 
+        /// Construct new cubic spline interpolation.
         /// </summary>
         /// <param name="xValues">List of x-values.</param>
         /// <param name="yValues">List of y-values.</param>

Numerics/Data/Interpolation/Support/Search.cs

@@ -762,7 +762,7 @@ public static int Hunt(double xValue, OrderedPairedData orderedPairedData, int s
             int XLO = start;
             int XHI;
             int XM;
-            var ASCND = default(bool);
+            var ASCND = orderedPairedData.OrderX == SortOrder.Ascending;
             int INC;
 
             // validate
@@ -816,7 +816,7 @@ public static int Hunt(double xValue, OrderedPairedData orderedPairedData, int s
             XLO = start;
 
             // Perform the hunt search algorithm
-            if (XLO <= 0 | XLO > N)
+            if (XLO <= 0 || XLO > N)
             {
                 // The input guess is not useful.
                 // Go immediately to bisection. 
@@ -849,7 +849,6 @@ public static int Hunt(double xValue, OrderedPairedData orderedPairedData, int s
                 else
                 {
                     // Hunt down
-                    XHI = XLO;
                     XHI = XLO - 1;
                     while (X < orderedPairedData[XLO].X == ASCND)
                     {
@@ -965,7 +964,7 @@ public static int Hunt(double xValue, IList<Ordinate> ordinateData, int start =
             XLO = start;
 
             // Perform the hunt search algorithm
-            if (XLO <= 0 | XLO > N)
+            if (XLO <= 0 || XLO > N)
             {
                 // The input guess is not useful.
                 // Go immediately to bisection. 
@@ -998,7 +997,6 @@ public static int Hunt(double xValue, IList<Ordinate> ordinateData, int start =
                 else
                 {
                     // Hunt down
-                    XHI = XLO;
                     XHI = XLO - 1;
                     while (X < ordinateData[XLO].X == ASCND)
                     {

Numerics/Data/Paired Data/Ordinate.cs

@@ -68,7 +68,7 @@ public Ordinate(double xValue, double yValue)
             X = xValue;
             Y = yValue;
             IsValid = true;
-            if (double.IsInfinity(X) | double.IsNaN(X) || double.IsInfinity(Y) | double.IsNaN(Y)) { IsValid = false; }
+            if (double.IsInfinity(X) || double.IsNaN(X) || double.IsInfinity(Y) || double.IsNaN(Y)) { IsValid = false; }
         }
 
 

Numerics/Data/Statistics/Autocorrelation.cs

@@ -215,6 +215,7 @@ public enum Type
             if (acf == null) return null;
 
             double den = acf[0, 1];
+            if (den == 0) return null;
             for (int i = 0; i < acf.GetLength(0); i++)
                 acf[i, 1] /= den;
             return acf;
@@ -239,6 +240,7 @@ public enum Type
             if (acf == null) return null;
 
             double den = acf[0, 1];
+            if (den == 0) return null;
             for (int i = 0; i < acf.GetLength(0); i++)
                 acf[i, 1] /= den;
             return acf;

Numerics/Data/Statistics/BoxCox.cs

@@ -133,6 +133,7 @@ public static double LogJacobian(IList<double> values, double lambda)
         /// <param name="lambda">The transformation exponent. Range -5 to +5.</param>
         public static double Transform(double value, double lambda)
         {
+            if (value <= 0) return double.NaN;
             if (Math.Abs(lambda) > 5d) return double.NaN;
             if (Math.Abs(lambda) < 1e-8) return Math.Log(value);
             return (Math.Pow(value, lambda) - 1.0d) / lambda;

Numerics/Data/Statistics/GoodnessOfFit.cs

@@ -834,7 +834,7 @@ public static double RSR(IList<double> observedValues, IList<double> modeledValu
             {
                 variance += Tools.Sqr(observedValues[i] - obsMean);
             }
-            double stdDev = Math.Sqrt(variance / n);
+            double stdDev = Math.Sqrt(variance / (n - 1));
 
             return rmse / stdDev;
         }

Numerics/Data/Statistics/Histogram.cs

@@ -189,6 +189,7 @@ public Histogram(IList<double> data)
             // Get the bin boundaries
             LowerBound = data.Min();
             UpperBound = data.Max();
+            if (UpperBound == LowerBound) UpperBound = LowerBound + 1.0;
             BinWidth = (UpperBound - LowerBound) / NumberOfBins;
             // Add bins
             double xl = LowerBound;
@@ -217,6 +218,7 @@ public Histogram(IList<double> data, int numberOfBins)
             NumberOfBins = numberOfBins;
             LowerBound = data.Min();
             UpperBound = data.Max();
+            if (UpperBound == LowerBound) UpperBound = LowerBound + 1.0;
             BinWidth = (UpperBound - LowerBound) / NumberOfBins;
             // Add bins
             double xl = LowerBound;
@@ -322,6 +324,7 @@ public double Median
                 int total = 0;
                 for (int i = 0; i < _bins.Count; i++)
                     total += _bins[i].Frequency;
+                if (total == 0) return double.NaN;
                 int halfTotal = (int)(total / 2d);
                 int m = 0;
                 int v = 0;
@@ -458,7 +461,7 @@ public int GetBinIndexOf(double value)
             SortBins();
             if (value < _bins.First().LowerBound || value > _bins.Last().UpperBound)
             {
-                throw new ArgumentException("value", "The value is not contained with the histogram limits.");
+                throw new ArgumentException("The value is not contained with the histogram limits.", nameof(value));
             }
             int idx = Search.Bisection(value, _binLimits);
             return idx < 0 ? 0 : idx >= _binLimits.Count ? _binLimits.Count - 1 : idx;

Numerics/Data/Statistics/HypothesisTests.cs

@@ -65,6 +65,7 @@ public class HypothesisTests
         /// <returns>Returns the 2-sided p-value of the test statistic.</returns>
         public static double OneSampleTtest(IList<double> sample, double populationMean = 0d)
         {
+            if (sample.Count < 2) throw new ArgumentException("Sample must have at least 2 observations.", nameof(sample));
             var meanVar = Statistics.MeanVariance(sample);
             int N = sample.Count;
             double se = Math.Sqrt(meanVar.Item2) / Math.Sqrt(N);
@@ -82,6 +83,7 @@ public static double OneSampleTtest(IList<double> sample, double populationMean
         /// <returns>Returns the 2-sided p-value of the test statistic.</returns>
         public static double EqualVarianceTtest(IList<double> sample1, IList<double> sample2)
         {
+            if (sample1.Count + sample2.Count < 3) throw new ArgumentException("Combined sample size must be at least 3.");
             var meanVar1 = Statistics.MeanVariance(sample1);
             var meanVar2 = Statistics.MeanVariance(sample2);
             int N1 = sample1.Count;
@@ -144,10 +146,12 @@ public static double PairedTtest(IList<double> sample1, IList<double> sample2)
         /// <returns>Returns the p-value of the test statistic.</returns>
         public static double Ftest(IList<double> sample1, IList<double> sample2)
         {
+            if (sample1.Count < 2 || sample2.Count < 2) throw new ArgumentException("Each sample must have at least 2 observations.");
             var meanVar1 = Statistics.MeanVariance(sample1);
             var meanVar2 = Statistics.MeanVariance(sample2);
             int n1 = sample1.Count, n2 = sample2.Count;
             double ave1 = meanVar1.Item1, var1 = meanVar1.Item2, ave2 = meanVar2.Item1, var2 = meanVar2.Item2, df1, df2, f, pVal;
+            if (var1 == 0 && var2 == 0) return 1.0;
             if (var1 > var2)
             {
                 f = var1 / var2;
@@ -177,6 +181,8 @@ public static double Ftest(IList<double> sample1, IList<double> sample2)
         /// <param name="pValue">The p-value of the test statistic.</param>
         public static void FtestModels(double sseRestricted, double sseFull, int dfRestricted, int dfFull, out double fStat, out double pValue)
         {
+            if (dfRestricted == dfFull) throw new ArgumentException("Restricted and full model degrees of freedom cannot be equal.");
+            if (dfFull <= 0) throw new ArgumentException("Full model degrees of freedom must be positive.", nameof(dfFull));
             fStat = ((sseRestricted - sseFull) / (dfRestricted - dfFull)) / (sseFull / dfFull);
             pValue = 2.0 * Beta.Incomplete(0.5 * dfFull, 0.5 * dfRestricted, dfFull / (dfFull + dfRestricted * fStat));
             if (pValue > 1.0) pValue = 2d - pValue;
@@ -329,7 +335,6 @@ public static double LinearTrendTest(IList<double> indices, IList<double> sample
             var yVals = new Vector(sample.ToArray());
             var lm = new LinearRegression(xVals, yVals, true);
             var tdist = new StudentT(lm.DegreesOfFreedom);
-            double d = Math.Abs(lm.Parameters[1] / lm.ParameterStandardErrors[1]);
             return (1 - tdist.CDF(Math.Abs(lm.Parameters[1] / lm.ParameterStandardErrors[1]))) * 2;
         }
 

Numerics/Data/Statistics/MultipleGrubbsBeckTest.cs

@@ -173,10 +173,6 @@ public static int Function(double[] X)
             return MGBTP;
         }
 
-        private static int _nIn;
-        private static int _rIn;
-        private static double _etaIn;
-
         /// <summary>
         /// Auxiliary routine used to compute p-values (GGCRITP) for a Generalized Grubbs-Beck Test.
         /// </summary>
@@ -187,16 +183,10 @@ private static double GGBCRITP(int N, int R, double ETA)
             {
                 return 0.5d;
             }
-            else
-            {
-                _nIn = N;
-                _rIn = R;
-                _etaIn = ETA;
-            }
 
-            // The original FORTRAN source code utilized a globally adaptive Gauss-Kronrod integration method. 
+            // The original FORTRAN source code utilized a globally adaptive Gauss-Kronrod integration method.
             // The number of low outliers computed by this method is consistent with the results from the FORTRAN code.
-            var sr = new AdaptiveGaussKronrod(FGGB, 1E-16, 1 - 1E-16);
+            var sr = new AdaptiveGaussKronrod((pzr) => FGGB(pzr, N, R, ETA), 1E-16, 1 - 1E-16);
             sr.MaxDepth = 25;
             sr.ReportFailure = false;
             sr.Integrate();
@@ -206,16 +196,13 @@ private static double GGBCRITP(int N, int R, double ETA)
         /// <summary>
         /// Auxiliary routine used in GGBCRITP
         /// </summary>
-        private static double FGGB(double PZR)
+        private static double FGGB(double PZR, int N, int R, double ETA)
         {
             double df, MuM, MuS2, VarM, VarS2, CovMS2;
             double EX1, EX2, EX3, EX4;
             double CovMS, VarS, alpha, beta;
             double MuMP, EtaP, H, Lambda, MuS, ncp, q, VarMP, PR, ZR, N2;
             double ANS;
-            int N = _nIn;
-            int R = _rIn;
-            double ETA = _etaIn;
 
             // Compute the value of the r-th smallest obs. based on its order statistic
             N2 = N - R;
@@ -261,6 +248,7 @@ public static void GrubbsBeckTest(IList<double> sample, out double XHi, out doub
         {
             // The following polynomial approximation proposed by Pilon et al. (1985)
             int n = sample.Count;
+            if (sample.Any(x => x <= 0)) throw new ArgumentException("All sample values must be positive for the Grubbs-Beck test.", nameof(sample));
             var logSample = new double[n];
             for (int i = 0; i < n; i++) logSample[i] = Math.Log(sample[i]);
             double mean = Statistics.Mean(logSample);