Minor updates, and working on documentation.

Author: HadenSmith

Numerics/Data/Paired Data/ProbabilityOrdinate.cs

@@ -28,6 +28,8 @@
 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */
 
+using Numerics.Data;
+using Numerics.Data.Statistics;
 using Numerics.Mathematics.Optimization;
 using Numerics.Sampling.MCMC;
 using Numerics.Utilities;
@@ -56,6 +58,51 @@ namespace Numerics.Distributions
     public class UncertaintyAnalysisResults
     {
 
+        /// <summary>
+        /// Construct an instance of the UncertaintyAnalysisResults class.
+        /// </summary>
+        public UncertaintyAnalysisResults() { }
+
+        /// <summary>
+        /// Constructs a new instance with computed uncertainty metrics.
+        /// </summary>
+        /// <param name="parentDistribution">The parent (mode/point estimate) distribution.</param>
+        /// <param name="sampledDistributions">Array of sampled distributions from posterior or bootstrap.</param>
+        /// <param name="probabilities">Array of non-exceedance probabilities for quantile estimation.</param>
+        /// <param name="alpha">The confidence level (default = 0.1 for 90% CI).</param>
+        /// <param name="minProbability">Minimum probability for mean curve computation (default = 0.001).</param>
+        /// <param name="maxProbability">Maximum probability for mean curve computation (default = 1 - 1e-9).</param>
+        /// <param name="recordParameterSets">If true, stores all parameter sets from sampled distributions.</param>
+        public UncertaintyAnalysisResults(UnivariateDistributionBase parentDistribution,
+                                          UnivariateDistributionBase[] sampledDistributions,
+                                          double[] probabilities,
+                                          double alpha = 0.1,
+                                          double minProbability = 0.001,
+                                          double maxProbability = 1 - 1e-9,
+                                          bool recordParameterSets = false)
+        {
+            if (parentDistribution == null)
+                throw new ArgumentNullException(nameof(parentDistribution));
+            if (sampledDistributions == null || sampledDistributions.Length == 0)
+                throw new ArgumentException("Sampled distributions cannot be null or empty.", nameof(sampledDistributions));
+            if (probabilities == null || probabilities.Length == 0)
+                throw new ArgumentException("Probabilities cannot be null or empty.", nameof(probabilities));
+
+            ProcessModeCurve(parentDistribution, probabilities);
+            ProcessConfidenceIntervals(sampledDistributions, probabilities, alpha);
+            ProcessMeanCurve(sampledDistributions, probabilities, minProbability, maxProbability);
+
+            if (recordParameterSets)
+                ProcessParameterSets(sampledDistributions);
+
+            // Set default values
+            AIC = double.NaN;
+            BIC = double.NaN;
+            DIC = double.NaN;
+            RMSE = double.NaN;
+            ERL = double.NaN;
+        }
+
         /// <summary>
         /// The parent probability distribution.
         /// </summary>
@@ -302,32 +349,192 @@ public static UncertaintyAnalysisResults FromXElement(XElement xElement)
         }
 
         /// <summary>
-        /// Returns uncertainty analysis results for the distribution estimated from MCMC.
+        /// Process and set the parent distribution and computed curve (mode, plug-in, point estimate).
         /// </summary>
-        /// <param name="results">The MCMC results.</param>
-        /// <param name="distribution">The parent distribution.</param>
-        /// <param name="probabilities">List of non-exceedance probabilities.</param>
+        /// <param name="parentDistribution">The parent distribution.</param>
+        /// <param name="probabilities">Array of non-exceedance probabilities.</param>
+        public void ProcessModeCurve(UnivariateDistributionBase parentDistribution, double[] probabilities)
+        {
+            if (parentDistribution == null)
+                throw new ArgumentNullException(nameof(parentDistribution));
+            if (probabilities == null || probabilities.Length == 0)
+                throw new ArgumentException("Probabilities cannot be null or empty.", nameof(probabilities));
+
+            ParentDistribution = parentDistribution;
+            ModeCurve = ParentDistribution.InverseCDF(probabilities);
+        }
+
+        /// <summary>
+        /// Process and set the confidence intervals from a list of sampled distributions.
+        /// </summary>
+        /// <param name="sampledDistributions">The list of sampled distributions to process.</param>
+        /// <param name="probabilities">Array of non-exceedance probabilities.</param>
         /// <param name="alpha">The confidence level; Default = 0.1, which will result in the 90% confidence intervals.</param>
-        public static UncertaintyAnalysisResults FromMCMCResults(MCMCResults results, UnivariateDistributionBase distribution, IList<double> probabilities, double alpha = 0.1)
+        public void ProcessConfidenceIntervals(UnivariateDistributionBase[] sampledDistributions, double[] probabilities, double alpha = 0.1)
+        {
+            if (sampledDistributions == null || sampledDistributions.Length == 0)
+                throw new ArgumentException("Sampled distributions cannot be null or empty.", nameof(sampledDistributions));
+            if (probabilities == null || probabilities.Length == 0)
+                throw new ArgumentException("Probabilities cannot be null or empty.", nameof(probabilities));
+            if (alpha <= 0 || alpha >= 1)
+                throw new ArgumentOutOfRangeException(nameof(alpha), "Alpha must be between 0 and 1.");
+
+            int B = sampledDistributions.Length;
+            int p = probabilities.Length;
+            double lowerCI = alpha / 2d;
+            double upperCI = 1d - alpha / 2d;
+            ConfidenceIntervals = new double[p, 2];
+
+            // Loop over probabilities and record percentiles
+            for (int i = 0; i < p; i++)
+            {
+                var XValues = new double[B];
+
+                // Compute quantiles in parallel
+                Parallel.For(0, B, idx => {
+                    XValues[idx] = sampledDistributions[idx]?.InverseCDF(probabilities[i]) ?? double.NaN;
+                });
+
+                // Filter valid values and sort
+                int validCount = 0;
+                for (int j = 0; j < B; j++)
+                {
+                    if (!double.IsNaN(XValues[j])) validCount++;
+                }
+
+                var validValues = new double[validCount];
+                int writeIdx = 0;
+                for (int j = 0; j < B; j++)
+                {
+                    if (!double.IsNaN(XValues[j]))
+                        validValues[writeIdx++] = XValues[j];
+                }
+
+                Array.Sort(validValues);
+
+                // Record percentiles for CIs
+                ConfidenceIntervals[i, 0] = Statistics.Percentile(validValues, lowerCI, true);
+                ConfidenceIntervals[i, 1] = Statistics.Percentile(validValues, upperCI, true);
+            }
+        }
+
+        /// <summary>
+        /// Computes the mean (predictive) curve by averaging CDFs across all sampled distributions.
+        /// Uses log-spaced quantiles for efficient computation across wide ranges.
+        /// </summary>
+        /// <param name="sampledDistributions">Array of sampled distributions from posterior or bootstrap.</param>
+        /// <param name="probabilities">Array of non-exceedance probabilities for interpolation.</param>
+        /// <param name="minProbability">Minimum probability for range determination (default = 0.001).</param>
+        /// <param name="maxProbability">Maximum probability for range determination (default = 1 - 1e-9).</param>
+        public void ProcessMeanCurve(UnivariateDistributionBase[] sampledDistributions, double[] probabilities, double minProbability = 0.001, double maxProbability = 1 - 1e-9)
         {
-            if (results == null|| results.Output == null ||results.Output.Count == 0) return null;
-            if (results.Output[0].Values.Length != distribution.NumberOfParameters) return null;
+            if (sampledDistributions == null || sampledDistributions.Length == 0)
+                throw new ArgumentException("Sampled distributions cannot be null or empty.", nameof(sampledDistributions));
+            if (probabilities == null || probabilities.Length == 0)
+                throw new ArgumentException("Probabilities cannot be null or empty.", nameof(probabilities));
+
+            int B = sampledDistributions.Length;
 
-            int realz = results.Output.Count;
-            var dists = new UnivariateDistributionBase[results.Output.Count];
+            // Compute min and max X values across all distributions
+            double minX = double.MaxValue;
+            double maxX = double.MinValue;
+            object lockObject = new object();
 
-            Parallel.For(0, realz, idx => 
+            Parallel.For(0, B, j =>
             {
-                var d = distribution.Clone();
-                d.SetParameters(results.Output[idx].Values);
-                dists[idx] = d;
+                if (sampledDistributions[j] != null)
+                {
+                    var innerMin = sampledDistributions[j].InverseCDF(minProbability);
+                    var innerMax = sampledDistributions[j].InverseCDF(maxProbability);
+
+                    lock (lockObject)
+                    {
+                        if (innerMin < minX) minX = innerMin;
+                        if (innerMax > maxX) maxX = innerMax;
+                    }
+                }
             });
 
-            // Set up dummy bootstrap analysis
-            var boot = new BootstrapAnalysis(distribution, ParameterEstimationMethod.MaximumLikelihood, 100, realz);
-            var uar = boot.Estimate(probabilities, alpha, dists);
-            uar.ParameterSets = null;
-            return uar;
+            // Create log-spaced quantiles for efficient coverage
+            double shift = minX <= 0 ? Math.Abs(minX) + 1d : 0;
+            double min = minX + shift;
+            double max = maxX + shift;
+            int order = (int)Math.Floor(Math.Log10(max) - Math.Log10(min));
+            int bins = Math.Max(200, Math.Min(1000, 100 * order));
+
+            var quantiles = new double[bins];
+            double delta = (Math.Log10(max) - Math.Log10(min)) / (bins - 1);
+
+            for (int i = 0; i < bins; i++)
+            {
+                double logX = Math.Log10(min) + i * delta;
+                quantiles[i] = Math.Pow(10, logX) - shift;
+            }
+
+            // Compute expected probability for each quantile
+            var expected = new double[bins];
+            for (int i = 0; i < bins; i++)
+            {
+                double total = 0d;
+                Parallel.For(0, B, () => 0d, (j, loop, sum) =>
+                {
+                    if (sampledDistributions[j] != null)
+                    {
+                        sum += sampledDistributions[j].CDF(quantiles[i]);
+                    }
+                    return sum;
+                }, z => Tools.ParallelAdd(ref total, z));
+                expected[i] = total / B;
+            }
+
+            // Build monotonic interpolation points
+            var yVals = new List<double> { quantiles[0] };
+            var xVals = new List<double> { expected[0] };
+            double minY = quantiles[0];
+            double maxY = quantiles[0];
+
+            for (int i = 1; i < bins; i++)
+            {
+                if (expected[i] > xVals[xVals.Count - 1])
+                {
+                    minY = Math.Min(minY, quantiles[i]);
+                    maxY = Math.Max(maxY, quantiles[i]);
+                    yVals.Add(quantiles[i]);
+                    xVals.Add(expected[i]);
+                }
+            }
+
+            // Determine if log transform is appropriate
+            bool useLogTransform = minY > 0 && (Math.Log10(maxY) - Math.Log10(minY)) > 1;
+
+            // Interpolate mean curve at requested probabilities
+            var linint = new Linear(xVals, yVals)
+            {
+                XTransform = Transform.NormalZ,
+                YTransform = useLogTransform ? Transform.Logarithmic : Transform.None
+            };
+            MeanCurve = linint.Interpolate(probabilities);
+        }
+
+        /// <summary>
+        /// Processes and stores the parameter sets from all sampled distributions.
+        /// </summary>
+        /// <param name="sampledDistributions">Array of sampled distributions to extract parameters from.</param>
+        public void ProcessParameterSets(UnivariateDistributionBase[] sampledDistributions)
+        {
+            if (sampledDistributions == null || sampledDistributions.Length == 0)
+                throw new ArgumentException("Sampled distributions cannot be null or empty.", nameof(sampledDistributions));
+
+            int B = sampledDistributions.Length;
+            ParameterSets = new ParameterSet[B];
+
+            Parallel.For(0, B, idx =>
+            {
+                if (sampledDistributions[idx] != null)
+                {
+                    ParameterSets[idx] = new ParameterSet(sampledDistributions[idx].GetParameters, double.NaN);
+                }
+            });
         }
 
     }

Numerics/Distributions/Univariate/Uncertainty Analysis/UncertaintyAnalysisResults.cs

@@ -8,14 +8,14 @@
     <DebugType>embedded</DebugType>
 
     <Description>Numerics is a free and open-source library for .NET developed by USACE-RMC, providing a comprehensive set of methods and algorithms for numerical computations and statistical analysis.</Description>
-    <PackageProjectUrl>https://github.com/USArmy-Corps-of-Engineers-RMC/Numerics</PackageProjectUrl>
-    <RepositoryUrl>https://github.com/USArmy-Corps-of-Engineers-RMC/Numerics</RepositoryUrl>
+    <PackageProjectUrl>https://github.com/USACE-RMC/Numerics</PackageProjectUrl>
+    <RepositoryUrl>https://github.com/USACE-RMC/Numerics</RepositoryUrl>
     <RepositoryType>git</RepositoryType>
     <Title>Numerics</Title>
     <GeneratePackageOnBuild>False</GeneratePackageOnBuild>
     <Company>USACE-RMC</Company>
     <PackageReadmeFile>README.md</PackageReadmeFile>
-    <PackageTags>interpolation; regression; statistics; machine learning; probability distributions; bootstrap analysis; distribution fitting; MCMC sampling; optimization;</PackageTags>
+    <PackageTags>interpolation; regression; statistics; machine learning; probability distributions; bootstrap analysis; distribution fitting; Bayesian MCMC sampling; optimization;</PackageTags>
     <PackageLicenseFile>LICENSE</PackageLicenseFile>
     <PackageRequireLicenseAcceptance>True</PackageRequireLicenseAcceptance>
     <PackageId>RMC.Numerics</PackageId>
@@ -29,6 +29,9 @@
   <PropertyGroup>
     <!-- 1591 is missing xml comment for a public member. 1587 is an xml comment not placed on a valid member (in our case, mostly license text at top of files) -->
     <NoWarn>CS1591;CS1587</NoWarn>
+    <Version>2.0.0</Version>
+    <PackageReleaseNotes>This major update delivers broad improvements across the entire Numerics library, with substantial enhancements to data analysis, statistical modeling, numerical methods, distributions, optimization, machine learning utilities, and overall performance and reliability. The update reflects over a year of incremental development, code modernization, and refinement.</PackageReleaseNotes>
+    <AssemblyVersion>2.0.0.0</AssemblyVersion>
   </PropertyGroup>
 
 

Numerics/Numerics.csproj

@@ -29,9 +29,12 @@
 */
 
 using Microsoft.VisualStudio.TestTools.UnitTesting;
+using Microsoft.VisualStudio.TestTools.UnitTesting.Logging;
 using Numerics.Data.Statistics;
 using Numerics.Distributions;
 using System;
+using System.Reflection;
+using static Microsoft.ApplicationInsights.MetricDimensionNames.TelemetryContext;
 
 namespace Data.Statistics
 {
@@ -260,12 +263,14 @@ public void Test_MAPE()
         /// Test that MAPE throws exception when observed values contain zero.
         /// </summary>
         [TestMethod]
-        [ExpectedException(typeof(ArgumentException))]
         public void Test_MAPE_ZeroObserved()
         {
             var observed = new double[] { 2.5, 0.0, 2.1, 1.4 };
             var modeled = new double[] { 3.0, -0.5, 2.0, 1.5 };
-            GoodnessOfFit.MAPE(observed, modeled);
+            var ex = Assert.Throws<Exception>(() =>
+            {
+                GoodnessOfFit.MAPE(observed, modeled);
+            });
         }
 
         /// <summary>
@@ -303,12 +308,14 @@ public void Test_sMAPE_PerfectFit()
         /// Test that sMAPE throws exception when both observed and modeled values are zero.
         /// </summary>
         [TestMethod]
-        [ExpectedException(typeof(ArgumentException))]
         public void Test_sMAPE_BothZero()
         {
             var observed = new double[] { 2.5, 0.0, 2.1, 1.4 };
             var modeled = new double[] { 3.0, 0.0, 2.0, 1.5 };
-            GoodnessOfFit.sMAPE(observed, modeled);
+            var ex = Assert.Throws<Exception>(() =>
+            {
+                GoodnessOfFit.sMAPE(observed, modeled);
+            });
         }
 
         /// <summary>
@@ -760,36 +767,42 @@ public void Test_ADTest()
         /// Test that methods throw appropriate exceptions when observed and modeled arrays have different lengths.
         /// </summary>
         [TestMethod]
-        [ExpectedException(typeof(ArgumentOutOfRangeException))]
         public void Test_UnequalArrayLengths_RMSE()
         {
             var observed = new double[] { 1.0, 2.0, 3.0 };
             var modeled = new double[] { 1.0, 2.0 };
-            GoodnessOfFit.RMSE(observed, modeled);
+            var ex = Assert.Throws<Exception>(() =>
+            {
+                GoodnessOfFit.RMSE(observed, modeled);
+            });
         }
 
         /// <summary>
         /// Test that methods throw appropriate exceptions when observed and modeled arrays have different lengths.
         /// </summary>
         [TestMethod]
-        [ExpectedException(typeof(ArgumentOutOfRangeException))]
         public void Test_UnequalArrayLengths_NSE()
         {
             var observed = new double[] { 1.0, 2.0, 3.0 };
             var modeled = new double[] { 1.0, 2.0 };
-            GoodnessOfFit.NashSutcliffeEfficiency(observed, modeled);
+            var ex = Assert.Throws<Exception>(() =>
+            {
+                GoodnessOfFit.NashSutcliffeEfficiency(observed, modeled);
+            });
         }
 
         /// <summary>
         /// Test that statistical tests throw appropriate exceptions with insufficient data.
         /// </summary>
         [TestMethod]
-        [ExpectedException(typeof(ArgumentOutOfRangeException))]
         public void Test_InsufficientData_KS()
         {
             var observed = new double[] { };
             var norm = new Normal(0, 1);
-            GoodnessOfFit.KolmogorovSmirnov(observed, norm);
+            var ex = Assert.Throws<Exception>(() =>
+            {
+                GoodnessOfFit.KolmogorovSmirnov(observed, norm);
+            });
         }
 
         #endregion