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
|
module Distribution
module Normal
module Ruby_
class << self
# Return a Proc object which returns a random number drawn
# from the normal distribution with mean +mean+ and
# standard deviation +sigma+, i.e. from N(+mean+,+sigma+^2).
#
# == Arguments
# * +mean+ - mean of the normal distribution
# * +sigma+ - standard deviation, a strictly positive number
# * +seed+ - seed, an integer value to set the initial state
#
# == Reference
# * http://www.taygeta.com/random/gaussian.html
#
def rng(mean = 0, sigma = 1, seed = nil)
seed = Random.new_seed if seed.nil?
r = Random.new(seed)
returned = 0
y1 = 0
y2 = 0
lambda do
if returned == 0
begin
x1 = 2.0 * r.rand - 1.0
x2 = 2.0 * r.rand - 1.0
w = x1 * x1 + x2 * x2
end while (w >= 1.0)
w = Math.sqrt((-2.0 * Math.log(w)) / w)
y1 = x1 * w
y2 = x2 * w
returned = 1
y1 * sigma + mean
else
returned = 0
y2 * sigma + mean
end
end
end
# random number within a gaussian distribution X ~ N(0,1)
def rngu
rng(0, 1, nil)
end
# Normal probability density function (pdf)
# With x=0 and sigma=1
def pdf(x)
(1.0 / SQ2PI) * Math.exp(-(x**2 / 2.0))
end
# Normal cumulative distribution function (cdf).
#
# Returns the integral of normal distribution
# over (-Infty, z].
#
def cdf(z)
return 0.0 if z < -12
return 1.0 if z > 12
return 0.5 if z == 0.0
if z > 0.0
e = true
else
e = false
z = -z
end
z = z.to_f
z2 = z * z
t = q = z * Math.exp(-0.5 * z2) / SQ2PI
3.step(199, 2) do |i|
prev = q
t *= z2 / i
q += t
return(e ? 0.5 + q : 0.5 - q) if q <= prev
end
e ? 1.0 : 0.0
end
# Return the inverse CDF or P-value of the corresponding integral
def quantile(qn)
b = [1.570796288, 0.03706987906, -0.8364353589e-3,
-0.2250947176e-3, 0.6841218299e-5, 0.5824238515e-5,
-0.104527497e-5, 0.8360937017e-7, -0.3231081277e-8,
0.3657763036e-10, 0.6936233982e-12]
if qn < 0.0 || 1.0 < qn
$stderr.printf("Error : qn <= 0 or qn >= 1 in pnorm()!\n")
return 0.0
end
qn == 0.5 and return 0.0
w1 = qn
w3 = -Math.log(4.0 * w1 * (1.0 - w1))
w1 = b[0]
1.upto 10 do |i|
w1 += b[i] * w3**i
end
qn > 0.5 and return Math.sqrt(w1 * w3)
-Math.sqrt(w1 * w3)
end
alias_method :p_value, :quantile
end
end
end
end
|
