Class: SimpleRandom

Inherits:
Object
  • Object
show all
Defined in:
lib/simple-random.rb

Instance Method Summary collapse

Constructor Details

#initialize ⇒ SimpleRandom

Returns a new instance of SimpleRandom.



2
3
4
5
 
# File 'lib/simple-random.rb', line 2

def initialize
  @m_w = 521288629
  @m_z = 362436069
end

Instance Method Details

#beta(a, b) ⇒ Object 



95
96
97
98
99
100
 
# File 'lib/simple-random.rb', line 95

def beta(a, b)
  raise "Alpha and beta parameters must be positive. Received a = #{a} and b = #{b}." unless a > 0 && b > 0
  u = gamma(a, 1)
  v = gamma(b, 1)
  u / (u + v)
end

#cauchy(median, scale) ⇒ Object 



108
109
110
111
112
 
# File 'lib/simple-random.rb', line 108

def cauchy(median, scale)
  raise 'Scale must be positive' if scale <= 0

  median + scale * Math.tan(Math::PI * (uniform - 0.5))
end

#chi_square(degrees_of_freedom) ⇒ Object 



87
88
89
 
# File 'lib/simple-random.rb', line 87

def chi_square(degrees_of_freedom)
  gamma(0.5 * degrees_of_freedom, 2.0)
end

#dirichlet(*parameters) ⇒ Object 



129
130
131
132
133
 
# File 'lib/simple-random.rb', line 129

def dirichlet(*parameters)
  sample = parameters.map { |a| gamma(a, 1) }
  sum = sample.inject(0.0) { |sum, g| sum + g }
  sample.map { |g| g / sum }
end

#exponential(mean = 1) ⇒ Object

Get exponential random sample with specified mean



41
42
43
44
 
# File 'lib/simple-random.rb', line 41

def exponential(mean = 1)
  raise 'Mean must be positive' if mean <= 0
  -1.0 * mean * Math.log(uniform)
end

#gamma(shape, scale) ⇒ Object

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.



62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
 
# File 'lib/simple-random.rb', line 62

def gamma(shape, scale)
  if shape >= 1.0
    d = shape - 1.0 / 3.0
    c = 1 / ((9 * d) ** 0.5)
    while true
      v = 0.0
      while v <= 0.0
        x = normal
        v = 1.0 + c * x
      end
      v = v ** 3
      u = uniform
      if u < (1.0 - 0.0331 * (x ** 4)) || Math.log(u) < (0.5 * (x ** 2) + d * (1.0 - v + Math.log(v)))
        return scale * d * v
      end
    end
  elsif shape <= 0.0
    raise 'Shape must be positive'
  else
    g = gamma(shape + 1.0, 1.0)
    w = uniform
    return scale * g * (w ** (1.0 / shape))
  end
end

#inverse_gamma(shape, scale) ⇒ Object 



91
92
93
 
# File 'lib/simple-random.rb', line 91

def inverse_gamma(shape, scale)
  1.0 / gamma(shape, 1.0 / scale)
end

#laplace(mean, scale) ⇒ Object 



120
121
122
123
 
# File 'lib/simple-random.rb', line 120

def laplace(mean, scale)
  u = uniform
  mean + Math.log(2) + ((u < 0.5 ? 1 : -1) * scale * Math.log(u < 0.5 ? u : 1 - u))
end

#log_normal(mu, sigma) ⇒ Object 



125
126
127
 
# File 'lib/simple-random.rb', line 125

def log_normal(mu, sigma)
  Math.exp(normal(mu, sigma))
end

#normal(mean = 0.0, standard_deviation = 1.0) ⇒ Object

Sample normal distribution with given mean and standard deviation



35
36
37
38
 
# File 'lib/simple-random.rb', line 35

def normal(mean = 0.0, standard_deviation = 1.0)
  raise 'Invalid standard deviation' if standard_deviation <= 0
  mean + standard_deviation * ((-2.0 * Math.log(uniform)) ** 0.5) * Math.sin(2.0 * Math::PI * uniform)
end

#set_seed(*args) ⇒ Object 



7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
 
# File 'lib/simple-random.rb', line 7

def set_seed(*args)
  if args.size > 1
    @m_w = args.first.to_i if args.first.to_i != 0
    @m_z = args.last.to_i if args.last.to_i != 0
  elsif args.first.is_a?(Numeric)
    @m_w = args.first.to_i if args.first.to_i != 0
  elsif args.first.is_a?(Time)
    x = (args.first.to_f * 1000000).to_i
    @m_w = x >> 16
    @m_z = x % 4294967296 # 2 ** 32
  else
    x = (Time.now.to_f * 1000000).to_i
    @m_w = x >> 16
    @m_z = x % 4294967296 # 2 ** 32
  end

  @m_w %= 4294967296
  @m_z %= 4294967296
end

#student_t(degrees_of_freedom) ⇒ Object 



114
115
116
117
118
 
# File 'lib/simple-random.rb', line 114

def student_t(degrees_of_freedom)
  raise 'Degrees of freedom must be positive' if degrees_of_freedom <= 0

  normal / ((chi_square(degrees_of_freedom) / degrees_of_freedom) ** 0.5)
end

#triangular(lower, mode, upper) ⇒ Object

Get triangular random sample with specified lower limit, mode, upper limit



47
48
49
50
51
52
53
54
55
56
57
 
# File 'lib/simple-random.rb', line 47

def triangular(lower, mode, upper)
  raise 'Upper limit must be larger than lower limit' if upper < lower
  raise 'Mode must lie between the upper and lower limits' if (mode < lower || mode > upper)
  f_c = (mode - lower) / (upper - lower)
  uniform_rand_num = uniform
  if uniform_rand_num < f_c
    lower + Math.sqrt(uniform_rand_num * (upper - lower) * (mode - lower))
  else
    upper - Math.sqrt((1 - uniform_rand_num) * (upper - lower) * (upper - mode))
  end
end

#uniform(lower = 0, upper = 1) ⇒ Object

Produce a uniform random sample from the open interval (lower, upper). The method will not return either end point.



29
30
31
32
 
# File 'lib/simple-random.rb', line 29

def uniform(lower = 0, upper = 1)
  raise 'Invalid range' if upper <= lower
  ((get_unsigned_int + 1) * (upper - lower) / 4294967296.0) + lower
end

#weibull(shape, scale) ⇒ Object 



102
103
104
105
106
 
# File 'lib/simple-random.rb', line 102

def weibull(shape, scale)
  raise 'Shape and scale must be positive' if shape <= 0.0 || scale <= 0.0

  scale * ((-Math.log(uniform)) ** (1.0 / shape))
end