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

def cdf(aa,bb,sigma, epsilon=0.0001, alpha=2.5, max_iterations=100) # :nodoc:
  raise "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 {|i|
    d[i]=0.0 if a[i].nil?
    e[i]=1.0 if b[i].nil?
  }
  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
      
      if a[i]!=nil
        d[i]=uPhi((a[i]-sumc).quo(c[i-1,i-1]))
      end
     # puts "sumc:#{sumc}"
      
      if 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=varsum+f[m]**2
    n+=1
    error=alpha*Math::sqrt((varsum.quo(n) - (intsum.quo(n))**2).quo(n))
    end while(error>epsilon and n<max_iterations)
  
    f=intsum.quo(n)
  #p intsum
  #puts "f:#{f}, n:#{n}, error:#{error}"
  f
end

.iPhi(pr) ⇒ Object 



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

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