ocean.math.random.NormalSource source code

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