source: trunk/sources/HeuristicLab.Random/3.3/GammaDistributedRandom.cs @ 14407

Last change on this file since 14407 was 14407, checked in by bburlacu, 3 years ago

#2703: Implement gamma distribution after "A Simple Method for Generating Gamma Variables" - Marsaglia & Tsang, ACM Transactions on Mathematical Software, Vol. 26, No. 3, September 2000, Pages 363–372.

File size: 4.4 KB
Line 
1#region License Information
2/* HeuristicLab
3 * Copyright (C) 2002-2016 Heuristic and Evolutionary Algorithms Laboratory (HEAL)
4 *
5 * This file is part of HeuristicLab.
6 *
7 * HeuristicLab is free software: you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation, either version 3 of the License, or
10 * (at your option) any later version.
11 *
12 * HeuristicLab is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with HeuristicLab. If not, see <http://www.gnu.org/licenses/>.
19 */
20#endregion
21
22using System;
23using HeuristicLab.Common;
24using HeuristicLab.Core;
25using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
26using HeuristicLab.Random;
27
28namespace HeuristicLab.ExpressionGenerator {
29  /// <summary>
30  /// Gamma distribution implemented after
31  /// "A Simple Method for Generating Gamma Variables" - Marsaglia & Tsang
32  /// ACM Transactions on Mathematical Software, Vol. 26, No. 3, September 2000, Pages 363–372.
33  /// </summary>
34  [Item("GammaDistributedRandom", "A pseudo random number generator for gamma distributed random numbers.")]
35  [StorableClass]
36  public sealed class GammaDistributedRandom : Item, IRandom {
37    [Storable]
38    private double shape;
39    public double Shape {
40      get { return shape; }
41      set { shape = value; }
42    }
43
44    [Storable]
45    private double rate;
46    public double Rate {
47      get { return rate; }
48      set { rate = value; }
49    }
50
51    [Storable]
52    private readonly IRandom random;
53
54    public GammaDistributedRandom() {
55      random = new MersenneTwister();
56    }
57
58    public GammaDistributedRandom(IRandom random, double shape, double rate) {
59      this.random = random;
60      this.shape = shape;
61      this.rate = rate;
62    }
63
64    [StorableConstructor]
65    private GammaDistributedRandom(bool deserializing) : base(deserializing) { }
66
67    private GammaDistributedRandom(GammaDistributedRandom original, Cloner cloner) : base(original, cloner) {
68    }
69
70    public override IDeepCloneable Clone(Cloner cloner) {
71      return new GammaDistributedRandom(this, cloner);
72    }
73
74    public void Reset() {
75      random.Reset();
76    }
77
78    public void Reset(int seed) {
79      random.Reset(seed);
80    }
81
82    public int Next() {
83      throw new NotImplementedException();
84    }
85
86    public int Next(int maxVal) {
87      throw new NotImplementedException();
88    }
89
90    public int Next(int minVal, int maxVal) {
91      throw new NotImplementedException();
92    }
93
94    public double NextDouble() {
95      return NextDouble(random, shape, rate);
96    }
97
98    /// <summary>
99    /// <para>Sample a value from a gamma distribution.</para>
100    /// <para>Implementation of "A Simple Method for Generating Gamma Variables" - Marsaglia & Tsang
101    /// ACM Transactions on Mathematical Software, Vol. 26, No. 3, September 2000, Pages 363–372.</para>
102    /// </summary>
103    /// <param name="uniformRandom">A uniformly-distributed random number generator.</param>
104    /// <param name="shape">The shape (k, α) of the Gamma distribution. Range: α ≥ 0.</param>
105    /// <param name="rate">The rate or inverse scale (β) of the Gamma distribution. Range: β ≥ 0.</param>
106    /// <returns>A sample from a Gamma distributed random variable.</returns>
107    public static double NextDouble(IRandom uniformRandom, double shape, double rate) {
108      if (double.IsPositiveInfinity(rate)) {
109        return shape;
110      }
111      var a = 1d;
112      if (shape < 1) {
113        a = Math.Pow(uniformRandom.NextDouble(), 1 / shape);
114        shape += 1;
115      }
116      var d = shape - 1d / 3d;
117      var c = 1 / Math.Sqrt(9 * d);
118
119      for (;;) {
120        double v, x;
121        do {
122          x = NormalDistributedRandom.NextDouble(uniformRandom, 0, 1);
123          v = 1 + c * x;
124        } while (v <= 0);
125
126        v = v * v * v;
127        x = x * x; // save a multiplication below
128        var u = uniformRandom.NextDouble();
129        if (u < 1 - 0.0331 * x * x || Math.Log(u) < 0.5 * x + d * (1 - v + Math.Log(v)))
130          return a * d * v / rate;
131      }
132    }
133  }
134}
Note: See TracBrowser for help on using the repository browser.