class RandomMultivariateNormal

# File lib/random_multivariate_normal.rb, line 5
def initialize(seed = Random.new_seed)
  @rnd = Random.new(seed)
end

# File lib/random_multivariate_normal.rb, line 9
def rand(mean, cov)
  rand_multivariate_norm(@rnd, mean, cov)
end

# File lib/random_multivariate_normal.rb, line 49
def array_as_zero_matrix(size)
  Array.new(size) { Array.new(size, 0) }
end

# 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

# 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

# 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

# File lib/random_multivariate_normal.rb, line 21
def rand_norm(rnd)
  box_muller(rnd.rand, rnd.rand)[0]
end