Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/sources/HeuristicLab.Random/3.3/NormalDistributedRandom.cs @ 4243

Last change on this file since 4243 was 3997, checked in by gkronber, 14 years ago

Minor improvements concerning efficiency of symbolic expression tree encoding data structures and operators. #1073

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