Module: Distribution::NormalMultivariate
- Defined in:
- lib/distribution/normalmultivariate.rb
Overview
Calculate cdf and inverse cdf for Multivariate Distribution.
Class Method Summary collapse
-
.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.
- .iPhi(pr) ⇒ Object
- .uPhi(x) ⇒ Object
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 |