Module: Distribution::NormalMultivariate

Defined in:
lib/distribution/normalmultivariate.rb

Overview

Calculate cdf and inverse cdf for Multivariate Distribution.

Class Method Summary collapse

Class Method Details

.cdf(aa, bb, sigma, epsilon = 0.0001, alpha = 2.5, max_iterations = 100) ⇒ Object

Returns multivariate cdf distribution

  • a is the array of lower values
  • b is the array of higher values
  • s is an symmetric positive definite covariance matrix


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
 
# File 'lib/distribution/normalmultivariate.rb', line 9

def cdf(aa, bb, sigma, epsilon = 0.0001, alpha = 2.5, max_iterations = 100) # :nodoc:
  fail "Doesn't work yet"
  a = [nil] + aa
  b = [nil] + bb
  m = aa.size
  sigma = sigma.to_gsl if sigma.respond_to? :to_gsl

  cc = GSL::Linalg::Cholesky.decomp(sigma)
  c = cc.lower
  intsum = 0
  varsum = 0
  n = 0
  d = Array.new(m + 1, nil)
  e = Array.new(m + 1, nil)
  f = Array.new(m + 1, nil)
  (1..m).each do|i|
    d[i] = 0.0 if a[i].nil?
    e[i] = 1.0 if b[i].nil?
  end
  d[1] = uPhi(a[1].quo(c[0, 0])) unless d[1] == 0
  e[1] = uPhi(b[1].quo(c[0, 0])) unless e[1] == 1
  f[1] = e[1] - d[1]

  error = 1000
  begin
      w = (m + 1).times.collect { |_i| rand * epsilon }
      y = []
      (2..m).each do |i|
        y[i - 1] = iPhi(d[i - 1] + w[i - 1] * (e[i - 1] - d[i - 1]))
        sumc = 0
        (1..(i - 1)).each do |j|
          sumc += c[i - 1, j - 1] * y[j]
        end

        d[i] = uPhi((a[i] - sumc).quo(c[i - 1, i - 1])) unless a[i].nil?
        # puts "sumc:#{sumc}"

        unless b[i].nil?
          # puts "e[#{i}] :#{c[i-1,i-1]}"
          e[i] = uPhi((b[i] - sumc).quo(c[i - 1, i - 1]))
        end
        f[i] = (e[i] - d[i]) * f[i - 1]
      end
      intsum += intsum + f[m]
      varsum += f[m]**2
      n += 1
      error = alpha * Math.sqrt((varsum.quo(n) - (intsum.quo(n))**2).quo(n))
    end while(error > epsilon && n < max_iterations)

  f = intsum.quo(n)
  # p intsum
  # puts "f:#{f}, n:#{n}, error:#{error}"
  f
end

.iPhi(pr) ⇒ Object 



64
65
66
 
# File 'lib/distribution/normalmultivariate.rb', line 64

def iPhi(pr)
  Distribution::Normal.p_value(pr)
end

.uPhi(x) ⇒ Object 



68
69
70
 
# File 'lib/distribution/normalmultivariate.rb', line 68

def uPhi(x)
  Distribution::Normal.cdf(x)
end