class RandomMultivariateNormal
Constants
- VERSION
Public Class Methods
new(seed = Random.new_seed)
click to toggle source
# File lib/random_multivariate_normal.rb, line 5 def initialize(seed = Random.new_seed) @rnd = Random.new(seed) end
Public Instance Methods
rand(mean, cov)
click to toggle source
# File lib/random_multivariate_normal.rb, line 9 def rand(mean, cov) rand_multivariate_norm(@rnd, mean, cov) end
Private Instance Methods
array_as_zero_matrix(size)
click to toggle source
# File lib/random_multivariate_normal.rb, line 49 def array_as_zero_matrix(size) Array.new(size) { Array.new(size, 0) } end
box_muller(x, y)
click to toggle source
# File lib/random_multivariate_normal.rb, line 25 def box_muller(x, y) z1 = Math.sqrt(-2 * Math.log(x)) * Math.cos(2 * Math::PI * y) z2 = Math.sqrt(-2 * Math.log(x)) * Math.sin(2 * Math::PI * y) [z1, z2] end
cholesky_foctor(mat)
click to toggle source
# File lib/random_multivariate_normal.rb, line 31 def cholesky_foctor(mat) raise ArgumentError.new('Not symmetric matrix') unless mat.symmetric? raise ArgumentError.new('Not positive definite matrix') unless mat.eigen.d.each(:diagonal).all?{|d| d > 0.0} l = array_as_zero_matrix(mat.row_size) mat.each_with_index do |e, row, col| lij = if row == col Math.sqrt(mat[col, col] - col.times.inject(0){|sum, j| sum + (l[col][j] ** 2)}) elsif row > col (1.0/l[col][col]) * (mat[row, col] - col.times.inject(0){|sum, j| sum + l[row][j] * l[col][j]}) else 0.0 end l[row][col] = lij end Matrix[*l] end
rand_multivariate_norm(rnd, mu, cov)
click to toggle source
# File lib/random_multivariate_normal.rb, line 15 def rand_multivariate_norm(rnd, mu, cov) z = Vector[*mu.size.times.map{rand_norm(rnd)}] l = cholesky_foctor(cov) l*z+mu end
rand_norm(rnd)
click to toggle source
# File lib/random_multivariate_normal.rb, line 21 def rand_norm(rnd) box_muller(rnd.rand, rnd.rand)[0] end