-
Notifications
You must be signed in to change notification settings - Fork 21
/
rng.jl
143 lines (127 loc) · 3.84 KB
/
rng.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
## return a uniform random sample from the interval (a, b)
function rand_uniform(a, b)
a + rand()*(b - a)
end
## return a random sample from a normal (Gaussian) distribution
function rand_normal(mean, stdev)
if stdev <= 0.0
error("standard deviation must be positive")
end
u1 = rand()
u2 = rand()
r = sqrt( -2.0*log(u1) )
theta = 2.0*pi*u2
mean + stdev*r*sin(theta)
end
## return a random sample from an exponential distribution
function rand_exponential(mean)
if mean <= 0.0
error("mean must be positive")
end
-mean*log(rand())
end
## return a random sample from a gamma distribution
function rand_gamma(shape, scale)
if shape <= 0.0
error("Shape parameter must be positive")
end
if scale <= 0.0
error("Scale parameter must be positive")
end
## Implementation based on "A Simple Method for Generating Gamma Variables"
## by George Marsaglia and Wai Wan Tsang.
## ACM Transactions on Mathematical Software
## Vol 26, No 3, September 2000, pages 363-372.
if shape >= 1.0
d = shape - 1.0/3.0
c = 1.0/sqrt(9.0*d)
while true
x = rand_normal(0, 1)
v = 1.0 + c*x
while v <= 0.0
x = rand_normal(0, 1)
v = 1.0 + c*x
end
v = v*v*v
u = rand()
xsq = x*x
if u < 1.0 -.0331*xsq*xsq || log(u) < 0.5*xsq + d*(1.0 - v + log(v))
return scale*d*v
end
end
else
g = rand_gamma(shape+1.0, 1.0)
w = rand()
return scale*g*pow(w, 1.0/shape)
end
end
## return a random sample from a chi square distribution
## with the specified degrees of freedom
function rand_chi_square(dof)
rand_gamma(0.5, 2.0*dof)
end
## return a random sample from an inverse gamma random variable
function rand_inverse_gamma(shape, scale)
## If X is gamma(shape, scale) then
## 1/Y is inverse gamma(shape, 1/scale)
1.0 / rand_gamma(shape, 1.0 / scale)
end
## return a sample from a Weibull distribution
function rand_weibull(shape, scale)
if shape <= 0.0
error("Shape parameter must be positive")
end
if scale <= 0.0
error("Scale parameter must be positive")
end
scale * pow(-log(rand()), 1.0 / shape)
end
## return a random sample from a Cauchy distribution
function rand_cauchy(median, scale)
if scale <= 0.0
error("Scale parameter must be positive")
end
p = rand()
median + scale*tan(pi*(p - 0.5))
end
## return a random sample from a Student t distribution
function rand_student_t(dof)
if dof <= 0
error("Degrees of freedom must be positive")
end
## See Seminumerical Algorithms by Knuth
y1 = rand_normal(0, 1)
y2 = rand_chi_square(dof)
y1 / sqrt(y2 / dof)
end
## return a random sample from a Laplace distribution
## The Laplace distribution is also known as the double exponential distribution.
function rand_laplace(mean, scale)
if scale <= 0.0
error("Scale parameter must be positive")
end
u = rand()
if u < 0.5
retval = mean + scale*log(2.0*u)
else
retval = mean - scale*log(2*(1-u))
end
retval
end
## return a random sample from a log-normal distribution
function rand_log_normal(mu, sigma)
return exp(rand_normal(mu, sigma))
end
## return a random sample from a beta distribution
function rand_beta(a, b)
if a <= 0 || b <= 0
error("Beta parameters must be positive")
end
## There are more efficient methods for generating beta samples.
## However such methods are a little more efficient and much more complicated.
## For an explanation of why the following method works, see
## http://www.johndcook.com/distribution_chart.html#gamma_beta
u = rand_gamma(a, 1.0)
v = rand_gamma(b, 1.0)
u / (u + v)
end