Dotnet/FastRng
Dotnet/FastRng
1
0
Fork 0
Code Issues 2 Pull Requests Actions Packages Projects Releases 5 Wiki Activity

Fixed shape fitter

Browse Source
This commit is contained in:
Thorsten Sommer 2020-10-30 21:48:53 +01:00
parent c007c6aef0
commit b84b488025
1 changed files with 25 additions and 21 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

46
FastRng/Double/ShapeFitter.cs
View File

@ -2,51 +2,55 @@ using System;
using System.Runtime.CompilerServices;
using System.Threading;
using System.Threading.Tasks;
using FastRng.Double.Distributions;
namespace FastRng.Double
{
/// <summary>
/// ShapeFitter is a rejection sampler, cf. https://en.wikipedia.org/wiki/Rejection_sampling
/// </summary>
public sealed class ShapeFitter
{
private readonly double[] probabilities;
private readonly double[] samples;
private readonly IRandom rng;
private readonly ushort sampleSize;
private readonly double threshold;
public ShapeFitter(Func<double, double> shapeFunction, IRandom rng, ushort sampleSize = 50, double threshold = 0.99)
private readonly double max;
private readonly double sampleSize;
private readonly IDistribution uniform = new Uniform();
public ShapeFitter(Func<double, double> shapeFunction, IRandom rng, ushort sampleSize = 50)
{
this.rng = rng;
this.threshold = threshold;
this.sampleSize = sampleSize;
this.samples = new double[sampleSize];
this.probabilities = new double[sampleSize];
var sampleStepSize = 1.0 / sampleSize;
var sampleStepSize = 1.0d / sampleSize;
var nextStep = 0.0 + sampleStepSize;
var maxValue = 0.0d;
for (var n = 0; n < sampleSize; n++)
{
this.probabilities[n] = shapeFunction(nextStep);
if (this.probabilities[n] > maxValue)
maxValue = this.probabilities[n];
nextStep += sampleStepSize;
}
this.max = maxValue;
}
public async ValueTask<double> NextNumber(CancellationToken token = default)
{
while (!token.IsCancellationRequested)
{
var nextNumber = await this.rng.GetUniform(token);
var nextBucket = (int)Math.Floor(nextNumber * this.sampleSize);
// var firstInBucket = this.samples[nextBucket] == 0;
this.samples[nextBucket] += this.probabilities[nextBucket];
// if (firstInBucket) // TODO: Could be an option (optional)
// return nextNumber;
var x = await this.rng.GetUniform(token);
var nextBucket = (int)Math.Floor(x * this.sampleSize);
var threshold = this.probabilities[nextBucket];
var y = await this.rng.NextNumber(0.0d, this.max, this.uniform, token);
if (this.samples[nextBucket] >= this.threshold)
{
this.samples[nextBucket] -= this.threshold;
return nextNumber;
}
if(y > threshold)
continue;
return x;
}
return double.NaN;