PageRenderTime 493ms CodeModel.GetById 18ms RepoModel.GetById 1ms app.codeStats 0ms

/Numerics-v2/Core/Statistics/Distributions/IDeviateGenerator.cs

#
C# | 318 lines | 180 code | 93 blank | 45 comment | 37 complexity | bb0e6af7e977b2526eec14e707480c76 MD5 | raw file
  1. using System;
  2. using System.Collections.Generic;
  3. using System.Linq;
  4. using System.Text;
  5. namespace Meta.Numerics.Statistics.Distributions {
  6. internal interface IDeviateGenerator {
  7. double GetNext (Random rng);
  8. }
  9. internal static class DeviateGeneratorFactory {
  10. public static IDeviateGenerator GetNormalGenerator () {
  11. return(new BoxMullerRejectionNormalDeviateGenerator());
  12. }
  13. public static IDeviateGenerator GetCauchyGenerator () {
  14. return(new CauchyGenerator());
  15. }
  16. public static IDeviateGenerator GetGammaGenerator (double alpha) {
  17. // choose a gamma generator depending on whether the shape parameter is less than one or greater than one
  18. if (alpha <= 1.0) {
  19. return(new AhrensDieterLowAlphaGammaGenerator(alpha));
  20. } else {
  21. return(new MarsagliaTsangGammaGenerator(GetNormalGenerator(), alpha));
  22. }
  23. }
  24. public static IDeviateGenerator GetBetaGenerator (double alpha, double beta) {
  25. return (new BetaFromGammaGenerator(GetGammaGenerator(alpha), GetGammaGenerator(beta)));
  26. }
  27. }
  28. // Many normal generators are described and compared in Thomas, et al., "Gaussian Random Number Generators", ACM Computing Surveys 39 (2007)
  29. // Most of our normal generators are written based on the description there.
  30. // Timing for generating 10^7 deviates in ms:
  31. // Inverse CDF 12600 ms
  32. // Box-Muller 1060 ms
  33. // Box-Muller Rejection 920 ms
  34. // Ratio of Uniforms 1210 ms
  35. // Leva's R of U 1130 ms
  36. // Our measurements indicate that simple Box-Muller with rejection is slightly faster than Leva's ratio of uniforms algorithm, and
  37. // the review article indicates that they are equally good, so we stick with Box-Muller for now. We have not yet attempted to implement
  38. // the Zigurraut algorithm.
  39. // One problem with Box-Muller is that the saving of the next value makes it stateful and therefore not thread-safe (at least not without adding
  40. // locking logic that would slow it down considerable). On the other hand the Random class itself isn't thread-safe, so this seems pretty minor.
  41. internal class BoxMullerRejectionNormalDeviateGenerator : IDeviateGenerator {
  42. private bool haveNextDeviate = false;
  43. private double nextDeviate;
  44. public double GetNext (Random rng) {
  45. if (haveNextDeviate) {
  46. haveNextDeviate = false;
  47. return (nextDeviate);
  48. } else {
  49. // just as in Box-Mueller we need a point in the unit disc
  50. // but in order to avoid trig functions, we find it by
  51. // picking a point in the circumscribing square, then
  52. // rejecting it if it lies outside the unit circle
  53. double x, y, r2;
  54. do {
  55. // pick a point in the square (-1,1)^2
  56. x = 2.0 * rng.NextDouble() - 1.0;
  57. y = 2.0 * rng.NextDouble() - 1.0;
  58. // determine if it lies in the unit disc
  59. r2 = x * x + y * y;
  60. } while ((r2 > 1.0) || (r2 == 0.0));
  61. double f = Math.Sqrt(-2.0 * Math.Log(r2) / r2);
  62. // store one deviate
  63. nextDeviate = f * x;
  64. haveNextDeviate = true;
  65. // return the other
  66. return (f * y);
  67. }
  68. }
  69. }
  70. #if PAST
  71. public class BoxMullerNormalGenerator : IDeviateGenerator {
  72. private bool haveNextDeviate = false;
  73. private double nextDeviate;
  74. public double GetNext (Random rng) {
  75. if (haveNextDeviate) {
  76. haveNextDeviate = false;
  77. return (nextDeviate);
  78. } else {
  79. // pick a point in the unit disc
  80. double u = rng.NextDouble();
  81. double t = 2.0 * Math.PI * rng.NextDouble();
  82. double a = Math.Sqrt(-2.0 * Math.Log(u));
  83. // store one deviate
  84. nextDeviate = a * Math.Sin(t);
  85. haveNextDeviate = true;
  86. // return the other
  87. return (a * Math.Cos(t));
  88. }
  89. }
  90. }
  91. public class RatioOfUniformsNormalGenerator : IDeviateGenerator {
  92. private static readonly double c0 = Math.Sqrt(2.0 / Math.E);
  93. private static readonly double c1 = 4.0 * Math.Pow(Math.E, 0.25);
  94. private static readonly double c2 = 4.0 * Math.Pow(Math.E, -1.35);
  95. public double GetNext (Random rng) {
  96. while (true) {
  97. // pic a point in (0,1) X (-\sqrt{2/e},\sqrt{2/2})
  98. double u = rng.NextDouble();
  99. double v = c0 * (2.0 * rng.NextDouble() - 1.0);
  100. // form ratio; this is our candidate deviate
  101. double x = v / u;
  102. // determine if the candidate is in acceptance or rejection region
  103. double x2 = x * x;
  104. if (x2 < 5.0 - c1 * u) {
  105. // inside squeeze, accept immediately
  106. return (x);
  107. } else if (x2 > c2 / u + 1.4) {
  108. // outside squeeze, reject immediately
  109. continue;
  110. }
  111. // between squeezes, do costly border evaluation
  112. if (v * v < -4.0 * u * u * Math.Log(u)) return (x);
  113. }
  114. }
  115. }
  116. public class LevaNormalGenerator : IDeviateGenerator {
  117. private static readonly double c = Math.Sqrt(2.0 / Math.E);
  118. private const double s = 0.449871, t = -0.386595;
  119. private const double a = 0.19600, b = 0.25472;
  120. private const double r1 = 0.27597, r2 = 0.27846;
  121. public double GetNext (Random rng) {
  122. while (true) {
  123. double u = rng.NextDouble();
  124. double v = c * (2.0 * rng.NextDouble() - 1.0);
  125. double x = u - s;
  126. double y = Math.Abs(v) - t;
  127. double q = x * x + y * (a * y - b * x);
  128. if (q < r1) {
  129. // inside squeeze, accept
  130. return (v / u);
  131. } else if (q > r2) {
  132. // outside squeeze, reject
  133. continue;
  134. }
  135. // evaluate exact border
  136. if (v * v < -4.0 * u * u * Math.Log(u)) return (v / u);
  137. }
  138. }
  139. }
  140. #endif
  141. // This Cauchy generator uses the same basic idea as Box-Muller with rejection. It is suggested in Numerical Recipies.
  142. // Our timing indicates that it is about 25% faster than simply evaluating the inverse tangent function of the inverse CDF.
  143. internal class CauchyGenerator : IDeviateGenerator {
  144. public double GetNext (Random rng) {
  145. // pick a random point within the unit semicircle using rejection
  146. double x, y;
  147. do {
  148. x = 2.0 * rng.NextDouble() - 1.0;
  149. y = rng.NextDouble();
  150. } while ((x * x + y * y > 1.0) || (y == 0.0));
  151. // its tangent is the tangent of a random angle and is thus Cauchy distributed
  152. return (x / y);
  153. }
  154. }
  155. // Using these rejection-based gamma generators is about six times faster than evaluating the inverse CDF.
  156. // The Ahrens-Dieter generator for \alpha < 1 is very well described by Best
  157. // Best also describes two improvements
  158. internal class AhrensDieterLowAlphaGammaGenerator : IDeviateGenerator {
  159. public AhrensDieterLowAlphaGammaGenerator (double alpha) {
  160. if (alpha > 1.0) throw new ArgumentOutOfRangeException("alpha");
  161. a = alpha;
  162. b = (Math.E + a) / Math.E;
  163. }
  164. double a, b;
  165. public double GetNext (Random rng) {
  166. while (true) {
  167. double P = b * rng.NextDouble();
  168. if (P < 1.0) {
  169. double x = Math.Pow(P, 1.0 / a);
  170. if (rng.NextDouble() > Math.Exp(-x)) continue;
  171. return (x);
  172. } else {
  173. double x = -Math.Log((b - P) / a);
  174. if (rng.NextDouble() > Math.Pow(x, a - 1.0)) continue;
  175. return (x);
  176. }
  177. }
  178. }
  179. }
  180. internal class MarsagliaTsangGammaGenerator : IDeviateGenerator {
  181. public MarsagliaTsangGammaGenerator (IDeviateGenerator normalGenerator, double alpha) {
  182. if (alpha < 1.0) throw new ArgumentOutOfRangeException("alpha");
  183. this.normalGenerator = normalGenerator;
  184. a = alpha;
  185. a1 = a - 1.0 / 3.0;
  186. a2 = 1.0 / Math.Sqrt(9.0 * a1);
  187. }
  188. IDeviateGenerator normalGenerator;
  189. double a, a1, a2;
  190. public double GetNext (Random rng) {
  191. double u, v, z;
  192. while (true) {
  193. // generate a candidate point
  194. do {
  195. z = normalGenerator.GetNext(rng);
  196. v = 1.0 + a2 * z;
  197. } while (v <= 0.0);
  198. v = v * v * v;
  199. u = rng.NextDouble();
  200. // check whether the point is outside the acceptance region
  201. // first via a simple interrior squeeze
  202. double z2 = z * z; double z4 = z2 * z2;
  203. if (u > 1.0 - 0.331 * z4) {
  204. // then, if necessary, via the exact rejection boundary
  205. if (Math.Log(u) > z2 / 2.0 + a1 * (1.0 - v + Math.Log(v))) continue;
  206. }
  207. return (a1 * v);
  208. }
  209. }
  210. }
  211. internal class BetaFromGammaGenerator : IDeviateGenerator {
  212. public BetaFromGammaGenerator (IDeviateGenerator alphaGenerator, IDeviateGenerator betaGenerator) {
  213. if (alphaGenerator == null) throw new ArgumentNullException("alphaGenerator");
  214. if (betaGenerator == null) throw new ArgumentNullException("betaGenerator");
  215. this.alphaGenerator = alphaGenerator;
  216. this.betaGenerator = betaGenerator;
  217. }
  218. private readonly IDeviateGenerator alphaGenerator, betaGenerator;
  219. public double GetNext (Random rng) {
  220. double x = alphaGenerator.GetNext(rng);
  221. double y = betaGenerator.GetNext(rng);
  222. return (x / (x + y));
  223. }
  224. }
  225. }