Implementing sampling for a distribution

(David Rutten) #1

I’m adding the shifted Gompertz distribution to my program in addition to the existing distributions already in MathNet. It’s all well and good but I don’t really know how to populate an array of doubles with random samples, as the `IContinuousDistribution.Samples` method demands.

Where do I even start with that? Clearly picking a random x and then testing whether another random value is below the pdf at that x would lead to 99% of samples being rejected. There must be a more efficient way.

``````/// <summary>
/// Evaluate the Shifted Gompertz Probability-Density-Function.
/// </summary>
/// <param name="b">Scale parameter.</param>
/// <param name="n">Shape parameter.</param>
/// <param name="x">X value.</param>
/// <returns>Probability at x.</returns>
private static double GompertzPdf(double b, double n, double x)
{
double ebx = Math.Exp(-b * x);
double p0 = b * ebx;
double p1 = Math.Exp(-n * ebx);
double p2 = 1 + n * (1 - ebx);
return p0 * p1 * p2;
}
/// <summary>
/// Evaluate the Shifted Gompertz Cumulative-Density-Function.
/// </summary>
/// <param name="b">Scale parameter.</param>
/// <param name="n">Shape parameter.</param>
/// <param name="x">X value.</param>
/// <returns>Cumulative probability at x.</returns>
private static double GompertzCdf(double b, double n, double x)
{
double ebx = Math.Exp(-b * x);
double p0 = 1 - ebx;
double p1 = Math.Exp(-n * ebx);
return p0 * p1;
}
``````

(Christoph Rüegg) #2

If we cannot leverage a transformation from related distributions, maybe we can compute the inverse CDF and use inverse transform sampling?

(David Rutten) #3

Possibly, if push comes to shove I can always approximate the cdf with a polyline and intersect it with horizontal rays. Not the fastest solution but probably good enough.

(David Rutten) #4

Yup, intersecting the cdf profile with horizontal rays at random values worked:

``````/// <summary>
/// Create an array of random samples in accordance with the Gompertz distribution.
/// </summary>
/// <param name="b">Scale parameter.</param>
/// <param name="n">Shape parameter.</param>
/// <param name="array">Array to populate.</param>
/// <param name="random">Random engine.</param>
private static void GompertzSamples(double b, double n, double[] array, Random random)
{
const int samples = 512;
var limit = GompertzFindLimit(b, n);
var factor = limit / samples;
var cdf = ArrayExtensionMethods.CreateArray(samples, i => GompertzCdf(b, n, i * factor));

for (int i = 0; i < array.Length; i++)
{
var r = random.NextDouble(cdf[0], cdf.Last());
if (cdf.InverseTransform(r, out int index, out var partial))
{
var x = (index + partial) * factor;
array[i] = x;
}
else
throw new InvalidOperationException("Inverse transform sampling of Gompertz cumulative density function failed.");
}
}
private static double GompertzFindLimit(double b, double n)
{
const double threshold0 = 0.9999999;
const double threshold1 = 0.99999999;
double upper = 8;
while (upper < 1000000)
if (GompertzCdf(b, n, upper) > threshold1)
break;
else
upper *= 2;

double lower = 0.0;
while (true)
{
double x = 0.5 * lower + 0.5 * upper;
double p = GompertzCdf(b, n, x);
if (p < threshold0)
lower = x;
else if (p > threshold1)
upper = x;
else
return x;
}
}
``````