Migrated gamma distribution to shape fitter

FastRng/Double/Distributions/Gamma.cs

@@ -1,69 +0,0 @@
-using System;
-using System.Threading;
-using System.Threading.Tasks;
-
-namespace FastRng.Double.Distributions
-{
-    public sealed class Gamma : IDistribution
-    {
-        private static readonly IDistribution NORMAL_DISTRIBUTION = new Normal();
-        
-        private double shape = 1.0;
-        
-        public IRandom Random { get; set; }
-
-        public double Shape
-        {
-            get => this.shape;
-            set
-            {
-                if(value <= 0.0)
-                    throw new ArgumentOutOfRangeException(message: "Shape must be greater than 0", null);
-                
-                this.shape = value;
-            }
-        }
-
-        public double Scale { get; set; } = 1.0;
-        
-        public async ValueTask<double> GetDistributedValue(CancellationToken token = default)
-        {
-            if (this.Random == null)
-                return double.NaN;
-            
-            // Implementation based on "A Simple Method for Generating Gamma Variables"
-            // by George Marsaglia and Wai Wan Tsang.  ACM Transactions on Mathematical Software
-            // Vol 26, No 3, September 2000, pages 363-372.
-            
-            if (this.Shape >= 1.0)
-            {
-                var d = this.Shape - 1.0 / 3.0;
-                var c = 1.0 / Math.Sqrt(9.0 * d);
-                while(true)
-                {
-                    double x, v;
-                    do
-                    {
-                        x = await this.Random.NextNumber(NORMAL_DISTRIBUTION, token);
-                        v = 1.0 + c * x;
-                    }
-                    while (v <= 0.0);
-                    
-                    v = v * v * v;
-                    var u = await this.Random.GetUniform(token);
-                    var xSquared = x * x;
-                    
-                    if (u < 1.0 - 0.0331 * xSquared * xSquared || Math.Log(u) < 0.5 * xSquared + d * (1.0 - v + Math.Log(v)))
-                        return this.Scale * d * v;
-                }
-            }
-            else
-            {
-                var dist = new Gamma{ Scale = 1, Shape = 1 + this.Shape};
-                var g = await this.Random.NextNumber(dist, token);
-                var w = await this.Random.GetUniform(token);
-                return this.Scale * g * Math.Pow(w, 1.0 / this.Shape);
-            }
-        }
-    }
-}

FastRng/Double/Distributions/GammaA5B15.cs

@@ -0,0 +1,45 @@
+using System;
+using System.Threading;
+using System.Threading.Tasks;
+
+namespace FastRng.Double.Distributions
+{
+    public sealed class GammaA5B15 : IDistribution
+    {
+        private const double ALPHA = 5.0;
+        private const double BETA = 15.0;
+        private const double CONSTANT = 0.341344210715475;
+
+        private static readonly double GAMMA_ALPHA;
+        private static readonly double BETA_TO_THE_ALPHA;
+        
+        private ShapeFitter fitter;
+        private IRandom random;
+
+        static GammaA5B15()
+        {
+            GAMMA_ALPHA = MathTools.Gamma(ALPHA);
+            BETA_TO_THE_ALPHA = Math.Pow(BETA, ALPHA);
+        }
+        
+        public IRandom Random
+        {
+            get => this.random;
+            set
+            {
+                this.random = value;
+                this.fitter = new ShapeFitter(GammaA5B15.ShapeFunction, this.random, 100);
+            }
+        }
+        
+        private static double ShapeFunction(double x) => CONSTANT * ((BETA_TO_THE_ALPHA * Math.Pow(x, ALPHA - 1.0d) * Math.Exp(-BETA * x)) / GAMMA_ALPHA);
+        
+        public async ValueTask<double> GetDistributedValue(CancellationToken token = default)
+        {
+            if (this.Random == null)
+                return double.NaN;
+            
+            return await this.fitter.NextNumber(token);
+        }
+    }
+}

FastRngTests/Double/Distributions/Gamma.cs

@@ -1,139 +0,0 @@
-using System;
-using System.Diagnostics.CodeAnalysis;
-using System.Linq;
-using System.Threading.Tasks;
-using FastRng.Double;
-using NUnit.Framework;
-
-namespace FastRngTests.Double.Distributions
-{
-    [ExcludeFromCodeCoverage]
-    public class Gamma
-    {
-        [Test]
-        [Category(TestCategories.COVER)]
-        [Category(TestCategories.NORMAL)]
-        public async Task TestGammaDistribution01()
-        {
-            const double SHAPE = 10.0;
-            const double SCALE = 2.0;
-            
-            const double MEAN = SHAPE * SCALE;
-            const double VARIANCE = SHAPE * SCALE * SCALE;
-            
-            var dist = new FastRng.Double.Distributions.Gamma{ Shape = SHAPE, Scale = SCALE };
-            var stats = new RunningStatistics();
-            var rng = new MultiThreadedRng();
-            
-            for (var n = 0; n < 100_000; n++)
-                stats.Push(await rng.NextNumber(dist));
-            
-            rng.StopProducer();
-            TestContext.WriteLine($"mean={MEAN} vs. {stats.Mean}");
-            TestContext.WriteLine($"variance={VARIANCE} vs {stats.Variance}");
-            
-            Assert.That(stats.Mean, Is.EqualTo(MEAN).Within(0.4), "Mean is out of range");
-            Assert.That(stats.Variance, Is.EqualTo(VARIANCE).Within(0.4), "Variance is out of range");
-        }
-        
-        [Test]
-        [Category(TestCategories.COVER)]
-        [Category(TestCategories.NORMAL)]
-        public async Task TestGammaDistribution02()
-        {
-            const double SHAPE = 0.5;
-            const double SCALE = 2.0;
-            
-            const double MEAN = SHAPE * SCALE;
-            const double VARIANCE = SHAPE * SCALE * SCALE;
-            
-            var dist = new FastRng.Double.Distributions.Gamma{ Shape = SHAPE, Scale = SCALE };
-            var stats = new RunningStatistics();
-            var rng = new MultiThreadedRng();
-            
-            for (var n = 0; n < 100_000; n++)
-                stats.Push(await rng.NextNumber(dist));
-            
-            rng.StopProducer();
-            TestContext.WriteLine($"mean={MEAN} vs. {stats.Mean}");
-            TestContext.WriteLine($"variance={VARIANCE} vs {stats.Variance}");
-            
-            Assert.That(stats.Mean, Is.EqualTo(MEAN).Within(0.1), "Mean is out of range");
-            Assert.That(stats.Variance, Is.EqualTo(VARIANCE).Within(0.1), "Variance is out of range");
-        }
-        
-        [Test]
-        [Category(TestCategories.COVER)]
-        [Category(TestCategories.NORMAL)]
-        public async Task TestGammaGeneratorWithRange03()
-        {
-            var rng = new MultiThreadedRng();
-            var dist = new FastRng.Double.Distributions.Gamma { Random = rng }; // Test default parameters
-            
-            var samples = new double[1_000];
-            for (var n = 0; n < samples.Length; n++)
-                samples[n] = await dist.GetDistributedValue();
-            
-            rng.StopProducer();
-            Assert.That(samples.Min(), Is.GreaterThanOrEqualTo(0.0), "Min is out of range");
-            Assert.That(samples.Max(), Is.LessThanOrEqualTo(1.0), "Max is out of range");
-        }
-        
-        [Test]
-        [Category(TestCategories.COVER)]
-        [Category(TestCategories.NORMAL)]
-        public async Task TestGammaGeneratorWithRange01()
-        {
-            var rng = new MultiThreadedRng();
-            var samples = new double[1_000];
-            for (var n = 0; n < samples.Length; n++)
-                samples[n] = await rng.NextNumber(-1.0, 1.0, new FastRng.Double.Distributions.Gamma());
-            
-            rng.StopProducer();
-            Assert.That(samples.Min(), Is.GreaterThanOrEqualTo(-1.0), "Min is out of range");
-            Assert.That(samples.Max(), Is.LessThanOrEqualTo(1.0), "Max is out of range");
-        }
-        
-        [Test]
-        [Category(TestCategories.COVER)]
-        [Category(TestCategories.NORMAL)]
-        public async Task TestGammaGeneratorWithRange02()
-        {
-            var rng = new MultiThreadedRng();
-            var samples = new double[1_000];
-            for (var n = 0; n < samples.Length; n++)
-                samples[n] = await rng.NextNumber(0.0, 1.0, new FastRng.Double.Distributions.Gamma());
-            
-            rng.StopProducer();
-            Assert.That(samples.Min(), Is.GreaterThanOrEqualTo(0.0), "Min is out of range");
-            Assert.That(samples.Max(), Is.LessThanOrEqualTo(1.0), "Max is out of range");
-        }
-
-        [Test]
-        [Category(TestCategories.COVER)]
-        [Category(TestCategories.NORMAL)]
-        public void ParameterTest01()
-        {
-            var dist = new FastRng.Double.Distributions.Gamma();
-            
-            Assert.DoesNotThrow(() => dist.Scale = -45);
-            Assert.DoesNotThrow(() => dist.Scale = 15);
-            Assert.DoesNotThrow(() => dist.Scale = 0);
-            
-            Assert.Throws<ArgumentOutOfRangeException>(() => dist.Shape = 0);
-            Assert.Throws<ArgumentOutOfRangeException>(() => dist.Shape = -78);
-            Assert.DoesNotThrow(() => dist.Shape = 0.0001);
-            Assert.DoesNotThrow(() => dist.Shape = 4);
-        }
-
-        [Test]
-        [Category(TestCategories.COVER)]
-        [Category(TestCategories.NORMAL)]
-        public async Task NoRandomNumberGenerator01()
-        {
-            var dist = new FastRng.Double.Distributions.Gamma();
-            Assert.DoesNotThrowAsync(async () => await dist.GetDistributedValue());
-            Assert.That(await dist.GetDistributedValue(), Is.NaN);
-        }
-    }
-}

FastRngTests/Double/Distributions/GammaA5B15.cs

@@ -0,0 +1,106 @@
+using System;
+using System.Diagnostics.CodeAnalysis;
+using System.Linq;
+using System.Threading.Tasks;
+using FastRng.Double;
+using NUnit.Framework;
+
+namespace FastRngTests.Double.Distributions
+{
+    [ExcludeFromCodeCoverage]
+    public class GammaA5B15
+    {
+        [Test]
+        [Category(TestCategories.COVER)]
+        [Category(TestCategories.NORMAL)]
+        public async Task TestGammaDistribution01()
+        {
+            var dist = new FastRng.Double.Distributions.GammaA5B15();
+            var fra = new FrequencyAnalysis();
+            var rng = new MultiThreadedRng();
+            
+            for (var n = 0; n < 100_000; n++)
+                fra.CountThis(await rng.NextNumber(dist));
+            
+            rng.StopProducer();
+            var result = fra.NormalizeAndPlotEvents(TestContext.WriteLine);
+
+            Assert.That(result[0], Is.EqualTo(0.0000929594237282).Within(0.0005));
+            Assert.That(result[1], Is.EqualTo(0.0012801746797876).Within(0.002));
+            Assert.That(result[2], Is.EqualTo(0.0055781488254349).Within(0.003));
+            
+            Assert.That(result[21], Is.EqualTo(0.9331608887752720).Within(0.09));
+            Assert.That(result[22], Is.EqualTo(0.9594734828891280).Within(0.09));
+            Assert.That(result[23], Is.EqualTo(0.9790895765535350).Within(0.09));
+            
+            Assert.That(result[50], Is.EqualTo(0.3478287795336570).Within(0.06));
+            
+            Assert.That(result[75], Is.EqualTo(0.0403399049422936).Within(0.009));
+            Assert.That(result[85], Is.EqualTo(0.0163628388658126).Within(0.009));
+            Assert.That(result[90], Is.EqualTo(0.0097147611446660).Within(0.005));
+            
+            Assert.That(result[97], Is.EqualTo(0.0041135143233153).Within(0.008));
+            Assert.That(result[98], Is.EqualTo(0.0036872732029996).Within(0.008));
+            Assert.That(result[99], Is.EqualTo(0.0033038503429554).Within(0.008));
+        }
+
+        [Test]
+        [Category(TestCategories.COVER)]
+        [Category(TestCategories.NORMAL)]
+        public async Task TestGammaGeneratorWithRange03()
+        {
+            var rng = new MultiThreadedRng();
+            var dist = new FastRng.Double.Distributions.GammaA5B15 { Random = rng }; // Test default parameters
+            
+            var samples = new double[1_000];
+            for (var n = 0; n < samples.Length; n++)
+                samples[n] = await dist.GetDistributedValue();
+            
+            rng.StopProducer();
+            Assert.That(samples.Min(), Is.GreaterThanOrEqualTo(0.0), "Min is out of range");
+            Assert.That(samples.Max(), Is.LessThanOrEqualTo(1.0), "Max is out of range");
+        }
+        
+        [Test]
+        [Category(TestCategories.COVER)]
+        [Category(TestCategories.NORMAL)]
+        public async Task TestGammaGeneratorWithRange01()
+        {
+            var rng = new MultiThreadedRng();
+            var dist = new FastRng.Double.Distributions.GammaA5B15();
+            var samples = new double[1_000];
+            for (var n = 0; n < samples.Length; n++)
+                samples[n] = await rng.NextNumber(-1.0, 1.0, dist);
+            
+            rng.StopProducer();
+            Assert.That(samples.Min(), Is.GreaterThanOrEqualTo(-1.0), "Min is out of range");
+            Assert.That(samples.Max(), Is.LessThanOrEqualTo(1.0), "Max is out of range");
+        }
+        
+        [Test]
+        [Category(TestCategories.COVER)]
+        [Category(TestCategories.NORMAL)]
+        public async Task TestGammaGeneratorWithRange02()
+        {
+            var rng = new MultiThreadedRng();
+            var dist = new FastRng.Double.Distributions.GammaA5B15();
+            var samples = new double[1_000];
+            for (var n = 0; n < samples.Length; n++)
+                samples[n] = await rng.NextNumber(0.0, 1.0, dist);
+            
+            rng.StopProducer();
+            Assert.That(samples.Min(), Is.GreaterThanOrEqualTo(0.0), "Min is out of range");
+            Assert.That(samples.Max(), Is.LessThanOrEqualTo(1.0), "Max is out of range");
+        }
+
+        [Test]
+        [Category(TestCategories.COVER)]
+        [Category(TestCategories.NORMAL)]
+        public async Task NoRandomNumberGenerator01()
+        {
+            var dist = new FastRng.Double.Distributions.GammaA5B15();
+            Assert.DoesNotThrowAsync(async () => await dist.GetDistributedValue());
+            Assert.That(await dist.GetDistributedValue(), Is.NaN);
+        }
+    }
+}