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