1  #region License Information


2  /* HeuristicLab


3  * Copyright (C) 20022016 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 


22  using System;


23  using HeuristicLab.Common;


24  using HeuristicLab.Core;


25  using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;


26  using HeuristicLab.Random;


27 


28  namespace 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 uniformlydistributed 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  }

