source: branches/HeuristicLab.Problems.GrammaticalOptimization/HeuristicLab.Common/Rand.cs @ 13273

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace HeuristicLab.Common {
  // code from Tom Minka's lightspeed toolbox (random.c)
  // original copyright: 
  // Written by Tom Minka
  // (c) Microsoft Corporation. All rights reserved.
  public static class Rand {

    /* Returns a sample from Normal(0,1) */

    public static double RandNormal(Random random) {
      double x, y, radius;
      /* Generate a random point inside the unit circle */
      do {
        x = 2 * random.NextDouble() - 1;
        y = 2 * random.NextDouble() - 1;
        radius = (x * x) + (y * y);
      } while ((radius >= 1.0) || (radius == 0.0));
      /* Box-Muller formula */
      radius = Math.Sqrt(-2 * Math.Log(radius) / radius);
      x *= radius;
      return x;
    }


    /* Returns a sample from Gamma(a, 1).
     * For Gamma(a,b), scale the result by b.
     */
    public static double GammaRand(Random random, double a) {
      /* Algorithm:
       * G. Marsaglia and W.W. Tsang, A simple method for generating gamma
       * variables, ACM Transactions on Mathematical Software, Vol. 26, No. 3,
       * Pages 363-372, September, 2000.
       * http://portal.acm.org/citation.cfm?id=358414
       */
      double boost, d, c, v;
      if (a < 1) {
        /* boost using Marsaglia's (1961) method: gam(a) = gam(a+1)*U^(1/a) */
        boost = Math.Exp(Math.Log(random.NextDouble()) / a);
        a++;
      } else boost = 1;
      d = a - 1.0 / 3; c = 1.0 / Math.Sqrt(9 * d);
      while (true) {
        double x, u;
        do {
          x = RandNormal(random);
          v = 1 + c * x;
        } while (v <= 0);
        v = v * v * v;
        x = x * x;
        u = random.NextDouble();
        if ((u < 1 - .0331 * x * x) ||
           (Math.Log(u) < 0.5 * x + d * (1 - v + Math.Log(v)))) break;
      }
      return (boost * d * v);
    }

    /* Returns a sample from Beta(a,b) */
    public static double BetaRand(Random random, double a, double b) {
      double g = GammaRand(random, a);
      return g / (g + GammaRand(random, b));
    }
  }
}