Add Bayesian probability comparison for perf results

Introduce Bayesian probability estimation using Student's t-distribution to compare session latencies. Add supporting statistical distribution classes and integrate results into DiffPrinter output. Extend unit tests to cover Bayesian comparison scenarios. Refactor variance calculation and update test dependencies to latest versions.

Author: ladeak

src/CHttp/Performance/Statitics/DiffPrinter.cs

@@ -37,7 +37,7 @@ public void Compare(PerformanceMeasurementResults session0, PerformanceMeasureme
 
         int lineLength = _console.WindowWidth;
         var scaleNormalize = (double)lineLength / (stats0.Durations.Length + stats0.Durations.Length);
-        string separator = new string('-', lineLength);
+        string separator = new('-', lineLength);
 
         // Histogram
         if (stats0.Durations.Length >= 100 && stats1.Durations.Length >= 100)
@@ -50,6 +50,9 @@ public void Compare(PerformanceMeasurementResults session0, PerformanceMeasureme
         PrintStatusCodes(stats0.StatusCodes, DiffStatusCodes(stats0, stats1));
         _console.WriteLine(separator);
 
+        if (PrintBayeisanComparison(session0, session1))
+            _console.WriteLine(separator);
+
         if (PrintWarnings(session0, session1))
             _console.WriteLine(separator);
     }
@@ -106,9 +109,7 @@ private static int[] DiffStatusCodes(Stats stats0, Stats stats1)
     {
         int[] diffStatusCodes = new int[stats0.StatusCodes.Length];
         for (int i = 0; i < diffStatusCodes.Length; i++)
-        {
             diffStatusCodes[i] = stats1.StatusCodes[i] - stats0.StatusCodes[i];
-        }
         return diffStatusCodes;
     }
 
@@ -172,4 +173,13 @@ private static long GetCountForBucket(Vector<double> bucketLimit, ref Span<long>
             currentCounter += input.Length;
         return currentCounter;
     }
+
+    private bool PrintBayeisanComparison(PerformanceMeasurementResults session0, PerformanceMeasurementResults session1)
+    {
+        double? percent = StatisticsCalculator.CalculateBayesianProbability(session0, session1) * 100;
+        if (percent.HasValue)
+            _console.WriteLine($"With {percent:00.0}% probability, the base session's true mean latency is lower than compared session's.");
+        return true;
+    }
+
 }

src/CHttp/Performance/Statitics/StatisticsCalculator.cs

@@ -1,4 +1,5 @@
-using System.Diagnostics.Metrics;
+using System.Diagnostics;
+using System.Diagnostics.Metrics;
 using System.Numerics;
 using CHttp.Data;
 using CHttp.Performance.Data;
@@ -45,7 +46,7 @@ public static Stats GetStats(PerformanceMeasurementResults session)
         Array.Sort(durations);
 
         var mean = totalTicks / (double)summaries.Count;
-        double stdDev = Math.Sqrt(CalcSquaredStdDev(durations, mean));
+        double stdDev = Math.Sqrt(CalculateVariance(durations, mean));
         double error = stdDev / Math.Sqrt(summaries.Count);
         double requestSec = (double)summaries.Count * TimeSpan.TicksPerSecond / (latestEnd - earliestStart);
         double throughput = session.TotalBytesRead / (mean / TimeSpan.TicksPerSecond);
@@ -72,7 +73,7 @@ public static Stats GetStats(PerformanceMeasurementResults session)
         return stats;
     }
 
-    private static double CalcSquaredStdDev(long[] durations, double mean)
+    private static double CalculateVariance(long[] durations, double mean)
     {
         var avg = new Vector<double>(mean);
         var input = durations.AsSpan();
@@ -86,7 +87,7 @@ private static double CalcSquaredStdDev(long[] durations, double mean)
                 var difference = Vector.Subtract(avg, vInput);
                 var squared = Vector.Multiply(difference, difference);
                 sum += Vector.Sum(squared);
-                input = input.Slice(vSize);
+                input = input[vSize..];
             }
         }
 
@@ -95,7 +96,32 @@ private static double CalcSquaredStdDev(long[] durations, double mean)
         {
             var difference = mean - input[0];
             sum += difference * difference;
-            input = input.Slice(1);
+            input = input[1..];
+        }
+
+        return sum / durations.Length;
+    }
+
+    private static double CalculateVariance(double[] durations, double mean)
+    {
+        var avg = new Vector<double>(mean);
+        var input = durations.AsSpan();
+        double sum = 0;
+        var vSize = Vector<double>.Count;
+        while (input.Length >= vSize)
+        {
+            var vInput = new Vector<double>(input);
+            var difference = Vector.Subtract(avg, vInput);
+            var squared = Vector.Multiply(difference, difference);
+            sum += Vector.Sum(squared);
+            input = input[vSize..];
+        }
+
+        // Remaining
+        for (int i = 0; i < input.Length; i++)
+        {
+            var difference = mean - input[i];
+            sum += difference * difference;
         }
 
         return sum / durations.Length;
@@ -148,4 +174,145 @@ public static (double DisplayValue, string Qualifier) Display(double value)
         }
         return (displayAverage, qualifier);
     }
+
+    internal static double? CalculateBayesianProbability(
+        PerformanceMeasurementResults session0,
+        PerformanceMeasurementResults session1,
+        Random? random = null)
+    {
+        if (session0.Summaries.Count < 100 || session0.Summaries.Count != session1.Summaries.Count)
+            return null;
+
+        // Step1: Calculate Mean and Variance based logarithmic durations
+        var (s0Mean, s0Variance) = LogMeanVariance(session0.Summaries);
+        var (s1Mean, s1Variance) = LogMeanVariance(session1.Summaries);
+
+        // The above values are not the "true" values, as these are based on the approximation of the samples.
+        // To calculate the "true" values use a Posterior Distributions (one for each session), in this case
+        // Student's t distribution ref: https://en.wikipedia.org/wiki/Student%27s_t-distribution
+        // and sample values of the this distribution. Compare the samples to derive the probability of one
+        // session being faster to the other session.
+
+        // For Student t distribution: number of requests, mean and variance needed.
+        random = random ?? Random.Shared;
+        var student0 = new StudentDistribution(session0.Summaries.Count - 1, s0Mean, Math.Sqrt(s0Variance / session0.Summaries.Count), random);
+        var student1 = new StudentDistribution(session1.Summaries.Count - 1, s1Mean, Math.Sqrt(s1Variance / session1.Summaries.Count), random);
+
+        // Simulate and compare
+        long session0Faster = 0;
+        double simulationCount = 8192;
+        for (int i = 0; i < simulationCount; i++)
+        {
+            if (student0.Sample() < student1.Sample())
+                session0Faster++;
+        }
+        return session0Faster / simulationCount;
+
+        static (double Mean, double Variance) LogMeanVariance(IReadOnlyCollection<Summary> summaries)
+        {
+            var durations = new double[summaries.Count];
+            double totalLogTicks = 0;
+
+            int current = 0;
+            foreach (var item in summaries)
+            {
+                var logTicks = Math.Log(item.Duration.Ticks);
+                durations[current++] = logTicks;
+                totalLogTicks += logTicks;
+            }
+            var mean = totalLogTicks / summaries.Count;
+            var variance = CalculateVariance(durations, mean);
+            return (mean, variance);
+        }
+    }
 }
+
+internal class StudentDistribution
+{
+    private readonly double _mean;
+    private readonly double _scale;
+    private readonly int _degreeOfFreedom;
+    private readonly StandardNormalDistribution _standardNormal;
+    private readonly ChiSquaredDistribution _chiSquared;
+
+    public StudentDistribution(int degreeOfFreedom, double mean, double scale, Random random)
+    {
+        _mean = mean;
+        _scale = scale;
+        _degreeOfFreedom = degreeOfFreedom;
+        _standardNormal = new(random);
+        _chiSquared = new(degreeOfFreedom, _standardNormal, random);
+    }
+
+    public double Sample()
+    {
+        // Student's t-distribution with v (degreeOfFreedom) degree of freedom can be defined as the distribution
+        // T = Z / Sqrt(V/degreeOfFreedom). 
+        // Z is a standard normal with expected value 0 and variance 1
+        // V is a chi-squared distribution
+        // degreeOfFreedom is the degree of freedom (sample count - 1)
+        var z = _standardNormal.Sample();
+        var v = _chiSquared.Sample();
+
+        var tStandard = z / Math.Sqrt(v / _degreeOfFreedom);
+
+        // Apply 'location' (mean) and 'scale' (variance).
+        return _mean + tStandard * _scale;
+    }
+}
+
+internal class StandardNormalDistribution(Random random)
+{
+    private readonly Random _random = random;
+    public double Sample()
+    {
+        // Ref: https://en.wikipedia.org/wiki/Normal_distribution
+        // The Box–Muller method uses two independent random numbers u1 and u2 distributed uniformly on (0,1).
+        // Then the two random variables X and Y
+        double u1 = 1.0 - _random.NextDouble();
+        double u2 = 1.0 - _random.NextDouble();
+        var sample = Math.Sqrt(-2.0 * Math.Log(u1)) * Math.Cos(2.0 * Math.PI * u2);
+        return sample;
+    }
+}
+
+internal class ChiSquaredDistribution(
+    int degreeOfFreedom,
+    StandardNormalDistribution standardNormal,
+    Random random)
+{
+    private readonly double _degreeOfFreedom = degreeOfFreedom;
+    private readonly Random _random = random;
+    public readonly StandardNormalDistribution _standardNormal = standardNormal;
+
+    public double Sample()
+    {
+        // Ref: https://en.wikipedia.org/wiki/Chi-squared_distribution
+        // Using the Gamma distribution method. x ~ Gamma(alpha = k/2, omega = 2)
+        // k in the above is (the degree of freedom, degreeOfFreedom variable).
+        return SampleGamma(_degreeOfFreedom / 2.0, 2);
+    }
+
+    private double SampleGamma(double alpha, double omega)
+    {
+        if (alpha < 1)
+            throw new InvalidOperationException("Unsupported Student-t distribution calculation");
+        // Ref: https://en.wikipedia.org/wiki/Gamma_distribution
+        // Marsaglia & Tsang GS algorithm, because this is faster
+        // to generating a sum of squared normal samples.
+        var d = alpha - 1 / 3.0;
+        var c = 1 / Math.Sqrt(9 * d);
+
+        while (true)
+        {
+            var z = _standardNormal.Sample();
+            var v = Math.Pow((1 + c * z), 3);
+            if (v <= 0)
+                continue;
+            var u = _random.NextDouble(); // Uniform distribution;
+
+            if (Math.Log(u) < 0.5 * z * z + d - d * v + d * Math.Log(v))
+                return d * v * omega;
+        }
+    }
+}

tests/CHttp.Api.Tests/CHttp.Api.Tests.csproj

@@ -19,7 +19,7 @@
   </ItemGroup>
 
   <ItemGroup>
-    <PackageReference Include="Microsoft.NET.Test.Sdk" Version="18.0.1" />
+    <PackageReference Include="Microsoft.NET.Test.Sdk" Version="18.6.0" />
     <PackageReference Include="xunit.runner.visualstudio" Version="3.1.5">
       <PrivateAssets>all</PrivateAssets>
       <IncludeAssets>runtime; build; native; contentfiles; analyzers; buildtransitive</IncludeAssets>

tests/CHttp.Tests/CHttp.Tests.csproj

@@ -9,9 +9,9 @@
   </PropertyGroup>
 
   <ItemGroup>
-    <PackageReference Include="Microsoft.AspNetCore.TestHost" Version="10.0.2" />
-    <PackageReference Include="Microsoft.Extensions.TimeProvider.Testing" Version="10.3.0" />
-    <PackageReference Include="Microsoft.NET.Test.Sdk" Version="18.0.1" />
+    <PackageReference Include="Microsoft.AspNetCore.TestHost" Version="11.0.0-preview.4.26230.115" />
+    <PackageReference Include="Microsoft.Extensions.TimeProvider.Testing" Version="10.6.0" />
+    <PackageReference Include="Microsoft.NET.Test.Sdk" Version="18.6.0" />
     <PackageReference Include="NSubstitute" Version="5.3.0" />
     <PackageReference Include="xunit.runner.visualstudio" Version="3.1.5">
       <PrivateAssets>all</PrivateAssets>