tango.math.random.NormalSource source code

1 /*******************************************************************************
2         copyright:      Copyright (c) 2008. Fawzi Mohamed
3         license:        BSD style: $(LICENSE)
4         version:        Initial release: July 2008
5         author:         Fawzi Mohamed
6 *******************************************************************************/
7 module tango.math.random.NormalSource;
8 private import Integer = tango.text.convert.Integer;
9 import tango.math.Math:exp,sqrt,log,PI;
10 import tango.math.ErrorFunction:erfc;
11 import tango.math.random.Ziggurat;
12 import tango.core.Traits: isRealType;
13 
14 /// class that returns gaussian (normal) distributed numbers (f=exp(-0.5*x*x)/sqrt(2*pi))
15 class NormalSource(RandG,T){
16     static assert(isRealType!(T),T.stringof~" not acceptable, only floating point variables supported");
17     /// probability distribution (non normalized, should be divided by sqrt(2*PI))
18     static real probDensityF(real x){ return exp(-0.5L*x*x); }
19     /// inverse probability distribution
20     static real invProbDensityF(real x){ return sqrt(-2.0L*log(x)); }
21     /// complement of the cumulative density distribution (integral x..infinity probDensityF)
22     static real cumProbDensityFCompl(real x){ return sqrt(0.5L*PI)*erfc(x/sqrt(2.0L));/*normalDistributionCompl(x);*/ }
23     /// tail for normal distribution
24     static T tailGenerator(RandG r, T dMin, int iNegative) 
25     { 
26         T x, y; 
27         do 
28         {
29             x = -log(r.uniform!(T)()) / dMin;
30             y = -log(r.uniform!(T)()); 
31         } while (y+y < x * x); 
32         return (iNegative ?(-x - dMin):(dMin + x)); 
33     }
34     alias Ziggurat!(RandG,T,probDensityF,tailGenerator,true) SourceT;
35     /// internal source of normal numbers
36     SourceT source;
37     /// initializes the probability distribution
38     this(RandG r){
39         source=SourceT.create!(invProbDensityF,cumProbDensityFCompl)(r,0xe.9dda4104d699791p-2L);
40     }
41     /// chainable call style initialization of variables (thorugh a call to randomize)
42     NormalSource opCall(U,S...)(ref U a,S args){
43         randomize(a,args);
44         return this;
45     }
46     /// returns a normal distribued number
47     final T getRandom(){
48         return source.getRandom();
49     }
50     /// returns a normal distribued number with the given sigma (standard deviation)
51     final T getRandom(T sigma){
52         return sigma*source.getRandom();
53     }
54     /// returns a normal distribued number with the given sigma (standard deviation) and mu (average)
55     final T getRandom(T sigma, T mu){
56         return mu+sigma*source.getRandom();
57     }
58     /// initializes a variable with normal distribued number and returns it
59     U randomize(U)(ref U a){
60         return source.randomize(a);
61     }
62     /// initializes a variable with normal distribued number with the given sigma and returns it
63     U randomize(U,V)(ref U a,V sigma){
64         return source.randomizeOp((T x){ return x*cast(T)sigma; },a);
65     }
66     /// initializes a variable with normal distribued numbers with the given sigma and mu and returns it
67     U randomize(U,V,S)(ref U el,V sigma, S mu){
68         return source.randomizeOp((T x){ return x*cast(T)sigma+cast(T)mu; },a);
69     }
70     /// initializes the variable with the result of mapping op on the random numbers (of type T)
71     U randomizeOp(U,S)(scope S delegate(T) op,ref U a){
72         return source.randomizeOp(op,a);
73     }
74     /// normal distribution with different default sigma and mu
75     /// f=exp(-x*x/(2*sigma^2))/(sqrt(2 pi)*sigma)
76     struct NormalDistribution{
77         T sigma,mu;
78         NormalSource source;
79         /// constructor
80         static NormalDistribution create(NormalSource source,T sigma,T mu){
81             NormalDistribution res;
82             res.source=source;
83             res.sigma=sigma;
84             res.mu=mu;
85             return res;
86         }
87         /// chainable call style initialization of variables (thorugh a call to randomize)
88         NormalDistribution opCall(U,S...)(ref U a,S args){
89             randomize(a,args);
90             return this;
91         }
92         /// returns a single number
93         T getRandom(){
94             return mu+sigma*source.getRandom();
95         }
96         /// initialize a
97         U randomize(U)(ref U a){
98             T op(T x){return mu+sigma*x; }
99             return source.randomizeOp(&op,a);
100         }
101         /// initialize a (let s and m have different types??)
102         U randomize(U,V)(ref U a,V s){
103             T op(T x){return mu+(cast(T)s)*x; }
104             return source.randomizeOp(&op,a);
105         }
106         /// initialize a (let s and m have different types??)
107         U randomize(U,V,S)(ref U a,V s, S m){
108             T op(T x){return (cast(T)m)+(cast(T)s)*x; }
109             return source.randomizeOp(&op,a);
110         }
111     }
112     /// returns a normal distribution with a non-default sigma/mu
113     NormalDistribution normalD(T sigma=cast(T)1,T mu=cast(T)0){
114         return NormalDistribution.create(this,sigma,mu);
115     }
116 }