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;
  }
}