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

Last change on this file since 11732 was 11732, checked in by gkronber, 5 years ago

#2283: refactoring and bug fixes

File size: 2.2 KB
Line 
1using System;
2using System.Collections.Generic;
3using System.Linq;
4using System.Text;
5using System.Threading.Tasks;
6
7namespace HeuristicLab.Common {
8  // code from Tom Minka's lightspeed toolbox (random.c)
9  // original copyright:
10  // Written by Tom Minka
11  // (c) Microsoft Corporation. All rights reserved.
12  public static class Rand {
13
14    /* Returns a sample from Normal(0,1) */
15
16    public static double RandNormal(Random random) {
17      double x, y, radius;
18      /* Generate a random point inside the unit circle */
19      do {
20        x = 2 * random.NextDouble() - 1;
21        y = 2 * random.NextDouble() - 1;
22        radius = (x * x) + (y * y);
23      } while ((radius >= 1.0) || (radius == 0.0));
24      /* Box-Muller formula */
25      radius = Math.Sqrt(-2 * Math.Log(radius) / radius);
26      x *= radius;
27      return x;
28    }
29
30
31    /* Returns a sample from Gamma(a, 1).
32     * For Gamma(a,b), scale the result by b.
33     */
34    public static double GammaRand(Random random, double a) {
35      /* Algorithm:
36       * G. Marsaglia and W.W. Tsang, A simple method for generating gamma
37       * variables, ACM Transactions on Mathematical Software, Vol. 26, No. 3,
38       * Pages 363-372, September, 2000.
39       * http://portal.acm.org/citation.cfm?id=358414
40       */
41      double boost, d, c, v;
42      if (a < 1) {
43        /* boost using Marsaglia's (1961) method: gam(a) = gam(a+1)*U^(1/a) */
44        boost = Math.Exp(Math.Log(random.NextDouble()) / a);
45        a++;
46      } else boost = 1;
47      d = a - 1.0 / 3; c = 1.0 / Math.Sqrt(9 * d);
48      while (true) {
49        double x, u;
50        do {
51          x = RandNormal(random);
52          v = 1 + c * x;
53        } while (v <= 0);
54        v = v * v * v;
55        x = x * x;
56        u = random.NextDouble();
57        if ((u < 1 - .0331 * x * x) ||
58           (Math.Log(u) < 0.5 * x + d * (1 - v + Math.Log(v)))) break;
59      }
60      return (boost * d * v);
61    }
62
63    /* Returns a sample from Beta(a,b) */
64    public static double BetaRand(Random random, double a, double b) {
65      double g = GammaRand(random, a);
66      return g / (g + GammaRand(random, b));
67    }
68  }
69}
Note: See TracBrowser for help on using the repository browser.