source: branches/3027-NormalDistribution/HeuristicLab.Random/3.3/NormalDistributedRandom.cs @ 17244

Last change on this file since 17244 was 17244, checked in by gkronber, 11 months ago

#3027: added class to generate normally distributed numbers using the polar method.

File size: 14.9 KB
Line 
1#region License Information
2/* HeuristicLab
3 * Copyright (C) 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 HEAL.Attic;
26
27namespace HeuristicLab.Random {
28
29  /// <summary>
30  /// Normally distributed random variable.
31  /// Uses the Ziggurat transform to efficiently generate normal distributed random numbers.
32  /// See "The Ziggurat Method for Generating Random Variables" (G. Marsaglia and W.W. Tsang 2000).
33  /// </summary>
34  [Item("NormalDistributedRandom", "A pseudo random number generator which uses the Ziggurat method to create normally distributed random numbers.")]
35  [StorableType("240F7050-C05B-4E85-82C3-2871FE2138B1")]
36  [Obsolete("Use NormalDistributedRandomPolar instead.")]
37  public sealed class NormalDistributedRandom : Item, IRandom {
38    [Storable]
39    private double mu;
40    /// <summary>
41    /// Gets or sets the value for µ.
42    /// </summary>
43    public double Mu {
44      get { return mu; }
45      set { mu = value; }
46    }
47
48    [Storable]
49    private double sigma;
50    /// <summary>
51    /// Gets or sets the value for sigma.
52    /// </summary>
53    public double Sigma {
54      get { return sigma; }
55      set { sigma = value; }
56    }
57
58    [Storable]
59    private IRandom uniform;
60
61    private static double[] w = new double[] {
62      1.7290404664e-09,
63      1.2680928529e-10,
64      1.6897518107e-10,
65      1.9862687883e-10,
66      2.2232431174e-10,
67      2.4244936614e-10,
68      2.6016130916e-10,
69      2.7611987696e-10,
70      2.9073962682e-10,
71      3.0429969655e-10,
72      3.1699795566e-10,
73      3.2898020419e-10,
74      3.4035738117e-10,
75      3.5121602848e-10,
76      3.6162509098e-10,
77      3.7164057942e-10,
78      3.8130856805e-10,
79      3.9066758162e-10,
80      3.9975012189e-10,
81      4.0858399997e-10,
82      4.1719308563e-10,
83      4.2559822333e-10,
84      4.3381759296e-10,
85      4.4186720949e-10,
86      4.4976131153e-10,
87      4.5751258337e-10,
88      4.6513240481e-10,
89      4.7263104541e-10,
90      4.8001774777e-10,
91      4.8730097735e-10,
92      4.9448850570e-10,
93      5.0158732723e-10,
94      5.0860404777e-10,
95      5.1554460700e-10,
96      5.2241466708e-10,
97      5.2921933502e-10,
98      5.3596349581e-10,
99      5.4265170135e-10,
100      5.4928817050e-10,
101      5.5587695558e-10,
102      5.6242188684e-10,
103      5.6892646150e-10,
104      5.7539412124e-10,
105      5.8182819673e-10,
106      5.8823168558e-10,
107      5.9460769641e-10,
108      6.0095900478e-10,
109      6.0728838625e-10,
110      6.1359850534e-10,
111      6.1989202660e-10,
112      6.2617133700e-10,
113      6.3243904558e-10,
114      6.3869737277e-10,
115      6.4494881657e-10,
116      6.5119559745e-10,
117      6.5744004685e-10,
118      6.6368432972e-10,
119      6.6993072201e-10,
120      6.7618144417e-10,
121      6.8243871665e-10,
122      6.8870464887e-10,
123      6.9498151678e-10,
124      7.0127148533e-10,
125      7.0757677495e-10,
126      7.1389966161e-10,
127      7.2024242126e-10,
128      7.2660727435e-10,
129      7.3299660786e-10,
130      7.3941280876e-10,
131      7.4585826404e-10,
132      7.5233547170e-10,
133      7.5884698525e-10,
134      7.6539541372e-10,
135      7.7198347714e-10,
136      7.7861395109e-10,
137      7.8528972214e-10,
138      7.9201378789e-10,
139      7.9878920145e-10,
140      8.0561923799e-10,
141      8.1250728368e-10,
142      8.1945689123e-10,
143      8.2647166888e-10,
144      8.3355555791e-10,
145      8.4071272166e-10,
146      8.4794732347e-10,
147      8.5526402627e-10,
148      8.6266754851e-10,
149      8.7016316375e-10,
150      8.7775620106e-10,
151      8.8545243360e-10,
152      8.9325818964e-10,
153      9.0117996399e-10,
154      9.0922497309e-10,
155      9.1740082198e-10,
156      9.2571583732e-10,
157      9.3417884539e-10,
158      9.4279972718e-10,
159      9.5158891877e-10,
160      9.6055785548e-10,
161      9.6971930486e-10,
162      9.7908692265e-10,
163      9.8867602993e-10,
164      9.9850361313e-10,
165      1.0085882129e-09,
166      1.0189509236e-09,
167      1.0296150599e-09,
168      1.0406069340e-09,
169      1.0519566329e-09,
170      1.0636980186e-09,
171      1.0758701707e-09,
172      1.0885182755e-09,
173      1.1016947354e-09,
174      1.1154610569e-09,
175      1.1298901814e-09,
176      1.1450695947e-09,
177      1.1611052120e-09,
178      1.1781275955e-09,
179      1.1962995039e-09,
180      1.2158286600e-09,
181      1.2369856250e-09,
182      1.2601323318e-09,
183      1.2857697129e-09,
184      1.3146201905e-09,
185      1.3477839955e-09,
186      1.3870635751e-09,
187      1.4357403044e-09,
188      1.5008658760e-09,
189      1.6030947680e-09
190    };
191    private static long[] k = new long[] {
192      1991057938,
193      0,
194      1611602771,
195      1826899878,
196      1918584482,
197      1969227037,
198      2001281515,
199      2023368125,
200      2039498179,
201      2051788381,
202      2061460127,
203      2069267110,
204      2075699398,
205      2081089314,
206      2085670119,
207      2089610331,
208      2093034710,
209      2096037586,
210      2098691595,
211      2101053571,
212      2103168620,
213      2105072996,
214      2106796166,
215      2108362327,
216      2109791536,
217      2111100552,
218      2112303493,
219      2113412330,
220      2114437283,
221      2115387130,
222      2116269447,
223      2117090813,
224      2117856962,
225      2118572919,
226      2119243101,
227      2119871411,
228      2120461303,
229      2121015852,
230      2121537798,
231      2122029592,
232      2122493434,
233      2122931299,
234      2123344971,
235      2123736059,
236      2124106020,
237      2124456175,
238      2124787725,
239      2125101763,
240      2125399283,
241      2125681194,
242      2125948325,
243      2126201433,
244      2126441213,
245      2126668298,
246      2126883268,
247      2127086657,
248      2127278949,
249      2127460589,
250      2127631985,
251      2127793506,
252      2127945490,
253      2128088244,
254      2128222044,
255      2128347141,
256      2128463758,
257      2128572095,
258      2128672327,
259      2128764606,
260      2128849065,
261      2128925811,
262      2128994934,
263      2129056501,
264      2129110560,
265      2129157136,
266      2129196237,
267      2129227847,
268      2129251929,
269      2129268426,
270      2129277255,
271      2129278312,
272      2129271467,
273      2129256561,
274      2129233410,
275      2129201800,
276      2129161480,
277      2129112170,
278      2129053545,
279      2128985244,
280      2128906855,
281      2128817916,
282      2128717911,
283      2128606255,
284      2128482298,
285      2128345305,
286      2128194452,
287      2128028813,
288      2127847342,
289      2127648860,
290      2127432031,
291      2127195339,
292      2126937058,
293      2126655214,
294      2126347546,
295      2126011445,
296      2125643893,
297      2125241376,
298      2124799783,
299      2124314271,
300      2123779094,
301      2123187386,
302      2122530867,
303      2121799464,
304      2120980787,
305      2120059418,
306      2119015917,
307      2117825402,
308      2116455471,
309      2114863093,
310      2112989789,
311      2110753906,
312      2108037662,
313      2104664315,
314      2100355223,
315      2094642347,
316      2086670106,
317      2074676188,
318      2054300022,
319      2010539237
320   };
321    private static double[] f = new double[] {
322      1.0000000000e+00,
323      9.6359968185e-01,
324      9.3628269434e-01,
325      9.1304361820e-01,
326      8.9228165150e-01,
327      8.7324303389e-01,
328      8.5550057888e-01,
329      8.3878362179e-01,
330      8.2290720940e-01,
331      8.0773830414e-01,
332      7.9317700863e-01,
333      7.7914607525e-01,
334      7.6558417082e-01,
335      7.5244158506e-01,
336      7.3967725039e-01,
337      7.2725689411e-01,
338      7.1515148878e-01,
339      7.0333611965e-01,
340      6.9178915024e-01,
341      6.8049186468e-01,
342      6.6942769289e-01,
343      6.5858197212e-01,
344      6.4794182777e-01,
345      6.3749545813e-01,
346      6.2723249197e-01,
347      6.1714339256e-01,
348      6.0721951723e-01,
349      5.9745317698e-01,
350      5.8783704042e-01,
351      5.7836467028e-01,
352      5.6902998686e-01,
353      5.5982738733e-01,
354      5.5075180531e-01,
355      5.4179835320e-01,
356      5.3296267986e-01,
357      5.2424055338e-01,
358      5.1562821865e-01,
359      5.0712203979e-01,
360      4.9871864915e-01,
361      4.9041482806e-01,
362      4.8220765591e-01,
363      4.7409430146e-01,
364      4.6607214212e-01,
365      4.5813870430e-01,
366      4.5029163361e-01,
367      4.4252872467e-01,
368      4.3484783173e-01,
369      4.2724698782e-01,
370      4.1972434521e-01,
371      4.1227802634e-01,
372      4.0490642190e-01,
373      3.9760786295e-01,
374      3.9038079977e-01,
375      3.8322380185e-01,
376      3.7613546848e-01,
377      3.6911445856e-01,
378      3.6215949059e-01,
379      3.5526937246e-01,
380      3.4844297171e-01,
381      3.4167915583e-01,
382      3.3497685194e-01,
383      3.2833510637e-01,
384      3.2175290585e-01,
385      3.1522938609e-01,
386      3.0876362324e-01,
387      3.0235484242e-01,
388      2.9600214958e-01,
389      2.8970485926e-01,
390      2.8346219659e-01,
391      2.7727350593e-01,
392      2.7113807201e-01,
393      2.6505529881e-01,
394      2.5902456045e-01,
395      2.5304529071e-01,
396      2.4711695313e-01,
397      2.4123899639e-01,
398      2.3541094363e-01,
399      2.2963231802e-01,
400      2.2390270233e-01,
401      2.1822164953e-01,
402      2.1258877218e-01,
403      2.0700371265e-01,
404      2.0146611333e-01,
405      1.9597564638e-01,
406      1.9053204358e-01,
407      1.8513499200e-01,
408      1.7978426814e-01,
409      1.7447963357e-01,
410      1.6922089458e-01,
411      1.6400785744e-01,
412      1.5884037316e-01,
413      1.5371830761e-01,
414      1.4864157140e-01,
415      1.4361007512e-01,
416      1.3862377405e-01,
417      1.3368265331e-01,
418      1.2878671288e-01,
419      1.2393598258e-01,
420      1.1913054436e-01,
421      1.1437050998e-01,
422      1.0965602100e-01,
423      1.0498725623e-01,
424      1.0036443919e-01,
425      9.5787845552e-02,
426      9.1257803142e-02,
427      8.6774669588e-02,
428      8.2338899374e-02,
429      7.7950984240e-02,
430      7.3611505330e-02,
431      6.9321118295e-02,
432      6.5080583096e-02,
433      6.0890771449e-02,
434      5.6752663106e-02,
435      5.2667401731e-02,
436      4.8636294901e-02,
437      4.4660862535e-02,
438      4.0742866695e-02,
439      3.6884389818e-02,
440      3.3087886870e-02,
441      2.9356317595e-02,
442      2.5693291798e-02,
443      2.2103304043e-02,
444      1.8592102453e-02,
445      1.5167297795e-02,
446      1.1839478277e-02,
447      8.6244847625e-03,
448      5.5489949882e-03,
449      2.6696291752e-03
450    };
451
452    [StorableConstructor]
453    private NormalDistributedRandom(StorableConstructorFlag _) : base(_) { }
454
455    /// <summary>
456    /// Initializes a new instance from an existing one (copy constructor).
457    /// </summary>
458    /// <param name="original">The original <see cref="NormalDistributedRandom"/> instance which is used to initialize the new instance.</param>
459    /// <param name="cloner">A <see cref="Cloner"/> which is used to track all already cloned objects in order to avoid cycles.</param>
460    private NormalDistributedRandom(NormalDistributedRandom original, Cloner cloner)
461      : base(original, cloner) {
462      uniform = cloner.Clone(original.uniform);
463      mu = original.mu;
464      sigma = original.sigma;
465    }
466
467    /// <summary>
468    /// Initializes a new instance of <see cref="NormalDistributedRandom"/> with µ = 0 and sigma = 1
469    /// and a new random number generator.
470    /// </summary>
471    public NormalDistributedRandom() {
472      this.mu = 0.0;
473      this.sigma = 1.0;
474      this.uniform = new MersenneTwister();
475    }
476
477    /// <summary>
478    /// Initializes a new instance of <see cref="NormalDistributedRandom"/> with the given parameters.
479    /// <note type="caution">The random number generator is not copied!</note>
480    /// </summary>   
481    /// <param name="uniformRandom">The random number generator.</param>
482    /// <param name="mu">The value for µ.</param>
483    /// <param name="sigma">The value for sigma.</param>
484    public NormalDistributedRandom(IRandom uniformRandom, double mu, double sigma) {
485      this.mu = mu;
486      this.sigma = sigma;
487      this.uniform = uniformRandom;
488    }
489
490    #region IRandom Members
491
492    /// <inheritdoc cref="IRandom.Reset()"/>
493    public void Reset() {
494      uniform.Reset();
495    }
496
497    /// <inheritdoc cref="IRandom.Reset(int)"/>
498    public void Reset(int seed) {
499      uniform.Reset(seed);
500    }
501
502    /// <summary>
503    /// This method is not implemented.
504    /// </summary>
505    public int Next() {
506      throw new NotImplementedException();
507    }
508
509    /// <summary>
510    /// This method is not implemented.
511    /// </summary>
512    public int Next(int maxVal) {
513      throw new NotImplementedException();
514    }
515
516    /// <summary>
517    /// This method is not implemented.
518    /// </summary>
519    public int Next(int minVal, int maxVal) {
520      throw new NotImplementedException();
521    }
522
523    /// <summary>
524    /// Generates a new double random number.
525    /// </summary>
526    /// <returns>A double random number.</returns>
527    public double NextDouble() {
528      return NormalDistributedRandom.NextDouble(uniform, mu, sigma);
529    }
530
531    #endregion
532
533    /// <summary>
534    /// Clones the current instance (deep clone).
535    /// </summary>
536    /// <remarks>Deep clone through <see cref="cloner.Clone"/> method of helper class.</remarks>
537    /// <returns>The cloned object as <see cref="NormalDistributedRandom"/>.</returns>
538    public override IDeepCloneable Clone(Cloner cloner) {
539      return new NormalDistributedRandom(this, cloner);
540    }
541
542    public static double NextDouble(IRandom uniformRandom, double mu, double sigma) {
543      double signFactor = uniformRandom.Next() % 2 == 0 ? 1.0 : -1.0;
544      return sigma * signFactor * NextPositiveDouble(uniformRandom) + mu;
545    }
546
547    private static double NextPositiveDouble(IRandom uniformRandom) {
548      int j = uniformRandom.Next();
549      int i = (j & 127);
550      if (Math.Abs(j) < k[i]) {
551        return j * w[i];
552      } else {
553        double r = 3.442620;
554        double x, y;
555        x = j * w[i];
556        if (i == 0) {
557          do {
558            x = -Math.Log(ScaledUniform(uniformRandom)) * 0.2904764;
559            y = -Math.Log(ScaledUniform(uniformRandom));
560          } while (y + y < x * x);
561          return (j > 0) ? r + x : -r - x;
562        }
563        if (f[i] + ScaledUniform(uniformRandom) * (f[i - 1] - f[i]) < Math.Exp(-0.5 * x * x)) {
564          return x;
565        } else {
566          // recurse
567          return (NextPositiveDouble(uniformRandom));
568        }
569      }
570    }
571
572    private static double ScaledUniform(IRandom uniformRandom) {
573      return 0.5 + uniformRandom.Next() * 0.2328306e-9;
574    }
575  }
576}
Note: See TracBrowser for help on using the repository browser.