1
2
3
4
5
6
7
8
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
65
66
67
68
69
70
71
72
73
|
module Distribution
# Calculate cdf and inverse cdf for Multivariate Distribution.
module NormalMultivariate
class << self
# 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
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
def iPhi(pr)
Distribution::Normal.p_value(pr)
end
def uPhi(x)
Distribution::Normal.cdf(x)
end
end
end
end
|
