From 74470f90bac33910fb25340eb71c993d0af552a7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=E5=9F=83=E5=8D=9A=E6=8B=89=E9=85=B1?= Date: Fri, 7 Jan 2022 14:27:40 +0800 Subject: [PATCH 1/5] RandomExtensions.NextBigIntegerSequence --- .../Random/RandomExtensionTests.cs | 18 +++++++++++ src/Numerics/Random/RandomExtensions.cs | 32 ++++++++++++++++++- 2 files changed, 49 insertions(+), 1 deletion(-) diff --git a/src/Numerics.Tests/Random/RandomExtensionTests.cs b/src/Numerics.Tests/Random/RandomExtensionTests.cs index 99c732c2..5bb84713 100644 --- a/src/Numerics.Tests/Random/RandomExtensionTests.cs +++ b/src/Numerics.Tests/Random/RandomExtensionTests.cs @@ -29,6 +29,8 @@ using MathNet.Numerics.Random; using NUnit.Framework; +using System.Numerics; +using System.Collections.Generic; namespace MathNet.Numerics.UnitTests.Random { @@ -48,6 +50,22 @@ namespace MathNet.Numerics.UnitTests.Random rnd.NextInt64(); } + /// + /// Can sample BigInteger + /// + [Test] + public void CanSampleBigInteger() + { + var rnd = new System.Random(0); + IEnumerator Sequence = rnd.NextBigIntegerSequence((BigInteger)long.MinValue * 3, (BigInteger)long.MaxValue * 3).GetEnumerator(); + Sequence.MoveNext(); + System.Console.WriteLine(Sequence.Current); + Sequence.MoveNext(); + System.Console.WriteLine(Sequence.Current); + Sequence.MoveNext(); + System.Console.WriteLine(Sequence.Current); + } + /// /// Can sample full range int32. /// diff --git a/src/Numerics/Random/RandomExtensions.cs b/src/Numerics/Random/RandomExtensions.cs index c351e3d8..0ba53982 100644 --- a/src/Numerics/Random/RandomExtensions.cs +++ b/src/Numerics/Random/RandomExtensions.cs @@ -29,6 +29,8 @@ using System; using System.Collections.Generic; +using System.Numerics; +using System.Linq; namespace MathNet.Numerics.Random { @@ -167,7 +169,7 @@ namespace MathNet.Numerics.Random } /// - /// Returns an infinite sequence of uniform random numbers greater than or equal to 0.0 and less than 1.0. + /// Returns an infinite sequence of uniform random 32-bit signed integers within the specified range. /// /// /// This extension is thread-safe if and only if called on an random number @@ -191,6 +193,34 @@ namespace MathNet.Numerics.Random } } + /// + /// Returns an infinite sequence of uniform random within the specified range. + /// + /// + /// This extension is thread-safe if and only if called on an random number + /// generator provided by Math.NET Numerics or derived from the RandomSource class. + /// + public static IEnumerable NextBigIntegerSequence(this System.Random rnd, BigInteger minInclusive, BigInteger maxExclusive) + { + BigInteger AbsoluteRange = maxExclusive - minInclusive; + int NumBytes = (int)Math.Ceiling(BigInteger.Log(AbsoluteRange, byte.MaxValue) * 2); + byte[] ByteSequence = Enumerable.Repeat(byte.MaxValue, NumBytes + 1).ToArray(); + ByteSequence[NumBytes] = 0; + BigInteger RandomNumber = new BigInteger(ByteSequence); + BigInteger ValidRange = RandomNumber - RandomNumber % AbsoluteRange; + while(true) + { + do + { + rnd.NextBytes(ByteSequence); + ByteSequence[NumBytes] = 0; + RandomNumber = new BigInteger(ByteSequence); + } + while (RandomNumber >= ValidRange); + yield return RandomNumber % AbsoluteRange + minInclusive; + } + } + /// /// Returns a nonnegative random number less than . /// From 2de37bbba159ad904259a140b0ad4fcf2b7005d3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=E5=9F=83=E5=8D=9A=E6=8B=89=E9=85=B1?= Date: Sat, 8 Jan 2022 15:03:00 +0800 Subject: [PATCH 2/5] Beta distribution may generate NaN when a and b are both very small --- .../DistributionTests/Continuous/BetaTests.cs | 8 ++++++++ src/Numerics/Distributions/Beta.cs | 2 +- 2 files changed, 9 insertions(+), 1 deletion(-) diff --git a/src/Numerics.Tests/DistributionTests/Continuous/BetaTests.cs b/src/Numerics.Tests/DistributionTests/Continuous/BetaTests.cs index 3351107c..fa28e6af 100644 --- a/src/Numerics.Tests/DistributionTests/Continuous/BetaTests.cs +++ b/src/Numerics.Tests/DistributionTests/Continuous/BetaTests.cs @@ -405,5 +405,13 @@ namespace MathNet.Numerics.UnitTests.DistributionTests.Continuous Assert.That(dist.InverseCumulativeDistribution(p), Is.EqualTo(x).Within(1e-6)); Assert.That(Beta.InvCDF(a, b, p), Is.EqualTo(x).Within(1e-6)); } + + [TestCase(0.001)] + public void ProbableNaNWhenABBothSmall(double ab) + { + Beta dist = new Beta(ab, ab); + for (byte i = 0; i < 100; ++i) + Assert.That(!double.IsNaN(dist.Sample()),"Generate NaN"); + } } } diff --git a/src/Numerics/Distributions/Beta.cs b/src/Numerics/Distributions/Beta.cs index 591317fe..8eec4c14 100644 --- a/src/Numerics/Distributions/Beta.cs +++ b/src/Numerics/Distributions/Beta.cs @@ -394,7 +394,7 @@ namespace MathNet.Numerics.Distributions { var x = Gamma.SampleUnchecked(rnd, a, 1.0); var y = Gamma.SampleUnchecked(rnd, b, 1.0); - return x/(x + y); + return x == 0 ? 0 : x / (x + y); } internal static void SamplesUnchecked(System.Random rnd, double[] values, double a, double b) From b9104ae50eb64a3c59755f8f4e27310da0e3cbaa Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=E5=9F=83=E5=8D=9A=E6=8B=89=E9=85=B1?= Date: Sat, 8 Jan 2022 15:11:01 +0800 Subject: [PATCH 3/5] Revert "RandomExtensions.NextBigIntegerSequence" This reverts commit 74470f90bac33910fb25340eb71c993d0af552a7. --- .../Random/RandomExtensionTests.cs | 18 ----------- src/Numerics/Random/RandomExtensions.cs | 32 +------------------ 2 files changed, 1 insertion(+), 49 deletions(-) diff --git a/src/Numerics.Tests/Random/RandomExtensionTests.cs b/src/Numerics.Tests/Random/RandomExtensionTests.cs index 5bb84713..99c732c2 100644 --- a/src/Numerics.Tests/Random/RandomExtensionTests.cs +++ b/src/Numerics.Tests/Random/RandomExtensionTests.cs @@ -29,8 +29,6 @@ using MathNet.Numerics.Random; using NUnit.Framework; -using System.Numerics; -using System.Collections.Generic; namespace MathNet.Numerics.UnitTests.Random { @@ -50,22 +48,6 @@ namespace MathNet.Numerics.UnitTests.Random rnd.NextInt64(); } - /// - /// Can sample BigInteger - /// - [Test] - public void CanSampleBigInteger() - { - var rnd = new System.Random(0); - IEnumerator Sequence = rnd.NextBigIntegerSequence((BigInteger)long.MinValue * 3, (BigInteger)long.MaxValue * 3).GetEnumerator(); - Sequence.MoveNext(); - System.Console.WriteLine(Sequence.Current); - Sequence.MoveNext(); - System.Console.WriteLine(Sequence.Current); - Sequence.MoveNext(); - System.Console.WriteLine(Sequence.Current); - } - /// /// Can sample full range int32. /// diff --git a/src/Numerics/Random/RandomExtensions.cs b/src/Numerics/Random/RandomExtensions.cs index 0ba53982..c351e3d8 100644 --- a/src/Numerics/Random/RandomExtensions.cs +++ b/src/Numerics/Random/RandomExtensions.cs @@ -29,8 +29,6 @@ using System; using System.Collections.Generic; -using System.Numerics; -using System.Linq; namespace MathNet.Numerics.Random { @@ -169,7 +167,7 @@ namespace MathNet.Numerics.Random } /// - /// Returns an infinite sequence of uniform random 32-bit signed integers within the specified range. + /// Returns an infinite sequence of uniform random numbers greater than or equal to 0.0 and less than 1.0. /// /// /// This extension is thread-safe if and only if called on an random number @@ -193,34 +191,6 @@ namespace MathNet.Numerics.Random } } - /// - /// Returns an infinite sequence of uniform random within the specified range. - /// - /// - /// This extension is thread-safe if and only if called on an random number - /// generator provided by Math.NET Numerics or derived from the RandomSource class. - /// - public static IEnumerable NextBigIntegerSequence(this System.Random rnd, BigInteger minInclusive, BigInteger maxExclusive) - { - BigInteger AbsoluteRange = maxExclusive - minInclusive; - int NumBytes = (int)Math.Ceiling(BigInteger.Log(AbsoluteRange, byte.MaxValue) * 2); - byte[] ByteSequence = Enumerable.Repeat(byte.MaxValue, NumBytes + 1).ToArray(); - ByteSequence[NumBytes] = 0; - BigInteger RandomNumber = new BigInteger(ByteSequence); - BigInteger ValidRange = RandomNumber - RandomNumber % AbsoluteRange; - while(true) - { - do - { - rnd.NextBytes(ByteSequence); - ByteSequence[NumBytes] = 0; - RandomNumber = new BigInteger(ByteSequence); - } - while (RandomNumber >= ValidRange); - yield return RandomNumber % AbsoluteRange + minInclusive; - } - } - /// /// Returns a nonnegative random number less than . /// From 1e0ca927f17162525f335797a631bf85f9b71f39 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=E5=9F=83=E5=8D=9A=E6=8B=89=E9=85=B1?= Date: Sat, 8 Jan 2022 15:20:59 +0800 Subject: [PATCH 4/5] Make Beta result equally possible to be 0 or 1 but not NaN when a and b are both very small --- src/Numerics/Distributions/Beta.cs | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/src/Numerics/Distributions/Beta.cs b/src/Numerics/Distributions/Beta.cs index 8eec4c14..c6e9f7c1 100644 --- a/src/Numerics/Distributions/Beta.cs +++ b/src/Numerics/Distributions/Beta.cs @@ -392,9 +392,13 @@ namespace MathNet.Numerics.Distributions /// a random number from the Beta distribution. internal static double SampleUnchecked(System.Random rnd, double a, double b) { - var x = Gamma.SampleUnchecked(rnd, a, 1.0); - var y = Gamma.SampleUnchecked(rnd, b, 1.0); - return x == 0 ? 0 : x / (x + y); + double x, y; + do + { + x = Gamma.SampleUnchecked(rnd, a, 1.0); + y = Gamma.SampleUnchecked(rnd, b, 1.0); + } while (x == 0 && y == 0); + return x / (x + y); } internal static void SamplesUnchecked(System.Random rnd, double[] values, double a, double b) From c578a385d440eae4eee731a9f1290696b961b963 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=E5=9F=83=E5=8D=9A=E6=8B=89=E9=85=B1?= Date: Sat, 8 Jan 2022 16:00:10 +0800 Subject: [PATCH 5/5] Handle a==b==0 case when it always generates NaN --- src/Numerics/Distributions/Beta.cs | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/src/Numerics/Distributions/Beta.cs b/src/Numerics/Distributions/Beta.cs index c6e9f7c1..67b07709 100644 --- a/src/Numerics/Distributions/Beta.cs +++ b/src/Numerics/Distributions/Beta.cs @@ -393,11 +393,20 @@ namespace MathNet.Numerics.Distributions internal static double SampleUnchecked(System.Random rnd, double a, double b) { double x, y; - do + if (a == b) { x = Gamma.SampleUnchecked(rnd, a, 1.0); y = Gamma.SampleUnchecked(rnd, b, 1.0); - } while (x == 0 && y == 0); + //When a==b (and possibly a==b==0), return value is equally possible to be 0 or 1 + if (x == 0 && y == 0) + return Bernoulli.Sample(0.5);//In particular, when a==b==0, Beta distribution degradates to Bernoulli distribution. + } + else + do + { + x = Gamma.SampleUnchecked(rnd, a, 1.0); + y = Gamma.SampleUnchecked(rnd, b, 1.0); + } while (x == 0 && y == 0);//When a!=b, return value is not equally possible to be 0 or 1. Regenerate. return x / (x + y); }