Class: Statsample::Bivariate::Polychoric::Processor

Inherits:
Object
  • Object
show all
Defined in:
lib/statsample/bivariate/polychoric/processor.rb

Overview

Provides statistics for a given combination of rho, alpha and beta and contingence table.

Constant Summary collapse

EPSILON = 
1e-10

Instance Attribute Summary collapse

Instance Method Summary collapse

Constructor Details

#initialize(alpha, beta, rho, matrix = nil) ⇒ Processor

Returns a new instance of Processor.



8
9
10
11
12
13
14
15
16
 
# File 'lib/statsample/bivariate/polychoric/processor.rb', line 8

def initialize(alpha,beta,rho,matrix=nil)
  @alpha=alpha
  @beta=beta
  @matrix=matrix
  @nr=@alpha.size+1
  @nc=@beta.size+1
  @rho=rho
  @pd=nil
end

Instance Attribute Details

#alpha ⇒ Object (readonly)

Returns the value of attribute alpha.



6
7
8
 
# File 'lib/statsample/bivariate/polychoric/processor.rb', line 6

def alpha
  @alpha
end

#beta ⇒ Object (readonly)

Returns the value of attribute beta.



6
7
8
 
# File 'lib/statsample/bivariate/polychoric/processor.rb', line 6

def beta
  @beta
end

#matrix ⇒ Object (readonly)

Returns the value of attribute matrix.



6
7
8
 
# File 'lib/statsample/bivariate/polychoric/processor.rb', line 6

def matrix
  @matrix
end

#rho ⇒ Object (readonly)

Returns the value of attribute rho.



6
7
8
 
# File 'lib/statsample/bivariate/polychoric/processor.rb', line 6

def rho
  @rho
end

Instance Method Details

#a(i) ⇒ Object 



37
38
39
40
 
# File 'lib/statsample/bivariate/polychoric/processor.rb', line 37

def a(i)
  raise "Index #{i} should be <= #{@nr-1}" if i>@nr-1
  i < 0 ? -100 : (i==@nr-1 ? 100 : alpha[i])
end

#b(j) ⇒ Object 



41
42
43
44
 
# File 'lib/statsample/bivariate/polychoric/processor.rb', line 41

def b(j)
  raise "Index #{j} should be <= #{@nc-1}" if j>@nc-1
  j < 0 ? -100 : (j==@nc-1 ? 100 : beta[j])
end

#bipdf(i, j) ⇒ Object 



18
19
20
 
# File 'lib/statsample/bivariate/polychoric/processor.rb', line 18

def bipdf(i,j)
   Distribution::BivariateNormal.pdf(a(i), b(j), rho)
end

#eq12(u, v) ⇒ Object 



46
47
48
 
# File 'lib/statsample/bivariate/polychoric/processor.rb', line 46

def eq12(u,v)
  Distribution::Normal.pdf(u)*Distribution::Normal.cdf((v-rho*u).quo( Math::sqrt(1-rho**2)))
end

#eq12b(u, v) ⇒ Object 



50
51
52
53
 
# File 'lib/statsample/bivariate/polychoric/processor.rb', line 50

def eq12b(u,v)
  Distribution::Normal.pdf(v) * Distribution::Normal.cdf((u-rho*v).quo( Math::sqrt(1-rho**2)))
  
end

#fd_loglike_a(k) ⇒ Object

First derivative for alpha_k Uses equation (6)



146
147
148
 
# File 'lib/statsample/bivariate/polychoric/processor.rb', line 146

def fd_loglike_a(k)
  fd_loglike_a_eq6(k)
end

#fd_loglike_a_eq13(k) ⇒ Object

Uses equation(13) from Olsson(1979)



169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
 
# File 'lib/statsample/bivariate/polychoric/processor.rb', line 169

def fd_loglike_a_eq13(k)
  rho=@rho
  if rho.abs>0.9999
    rho= (rho>0) ? 0.9999 : -0.9999
  end
  total=0
  a_k=a(k)
  @nc.times do |j|
    #puts "j: #{j}"
    #puts "b #{j} : #{b.call(j)}"
    #puts "b #{j-1} : #{b.call(j-1)}"
  
    e_1=@matrix[k,j].quo(pd[k][j]+EPSILON) - @matrix[k+1,j].quo(pd[k+1][j]+EPSILON)
    e_2=Distribution::Normal.pdf(a_k)
    e_3=Distribution::Normal.cdf((b(j)-rho*a_k).quo(Math::sqrt(1-rho**2))) - Distribution::Normal.cdf((b(j-1)-rho*a_k).quo(Math::sqrt(1-rho**2)))
    #puts "val #{j}: #{e_1} | #{e_2} | #{e_3}"
    total+= e_1*e_2*e_3
  end
  total
end

#fd_loglike_a_eq6(k) ⇒ Object

Uses equation (6) from Olsson(1979)



153
154
155
156
157
158
159
160
161
162
163
164
165
 
# File 'lib/statsample/bivariate/polychoric/processor.rb', line 153

def fd_loglike_a_eq6(k)
  rho=@rho
  if rho.abs>0.9999
    rho= (rho>0) ? 0.9999 : -0.9999
  end
  total=0
  @nr.times do |i|
    @nc.times  do |j|
      total+=@matrix[i,j].quo(pd[i][j]+EPSILON) * fd_loglike_cell_a(i,j,k)
    end
  end
  total
end

#fd_loglike_b(m) ⇒ Object

First derivative for beta_m. Uses equation 6 (Olsson,1979)



206
207
208
 
# File 'lib/statsample/bivariate/polychoric/processor.rb', line 206

def fd_loglike_b(m)
  fd_loglike_b_eq14(m)
end

#fd_loglike_b_eq14(m) ⇒ Object

First derivative for beta_m Uses equation(14) from Olsson(1979)



211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
 
# File 'lib/statsample/bivariate/polychoric/processor.rb', line 211

def fd_loglike_b_eq14(m)
  rho=@rho
  if rho.abs>0.9999
    rho= (rho>0) ? 0.9999 : -0.9999
  end
  total=0
  b_m=b(m)
  @nr.times do |i|
    e_1=@matrix[i,m].quo(pd[i][m]+EPSILON) - @matrix[i,m+1].quo(pd[i][m+1]+EPSILON)
    e_2=Distribution::Normal.pdf(b_m)
    e_3=Distribution::Normal.cdf((a(i)-rho*b_m).quo(Math::sqrt(1-rho**2))) - Distribution::Normal.cdf((a(i-1)-rho*b_m).quo(Math::sqrt(1-rho**2)))
    #puts "val #{j}: #{e_1} | #{e_2} | #{e_3}"
    
    total+= e_1*e_2*e_3
  end
  total
end

#fd_loglike_b_eq6(m) ⇒ Object

First derivative for b Uses equation 6 (Olsson, 1979)



191
192
193
194
195
196
197
198
199
200
201
202
203
 
# File 'lib/statsample/bivariate/polychoric/processor.rb', line 191

def fd_loglike_b_eq6(m)
  rho=@rho
  if rho.abs>0.9999
    rho= (rho>0) ? 0.9999 : -0.9999
  end
  total=0
  @nr.times do |i|
    @nc.times  do |j|
      total+=@matrix[i,j].quo(pd[i][j]+EPSILON) * fd_loglike_cell_b(i,j,m)
    end
  end
  total
end

#fd_loglike_cell_a(i, j, k) ⇒ Object

Equation(10) from Olsson(1979)



59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
 
# File 'lib/statsample/bivariate/polychoric/processor.rb', line 59

def fd_loglike_cell_a(i, j, k)
=begin
  if k==i
    Distribution::NormalBivariate.pd_cdf_x(a(k),b(j), rho) - Distribution::NormalBivariate.pd_cdf_x(a(k),b(j-1),rho)
  elsif k==(i-1)
    -Distribution::NormalBivariate.pd_cdf_x(a(k),b(j),rho) + Distribution::NormalBivariate.pd_cdf_x(a(k),b(j-1),rho)
  else
    0
  end
=end          
  if k==i
    eq12(a(k),b(j))-eq12(a(k), b(j-1))
  elsif k==(i-1)
    -eq12(a(k),b(j))+eq12(a(k), b(j-1))
  else
    0
  end          
end

#fd_loglike_cell_b(i, j, m) ⇒ Object 



78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
 
# File 'lib/statsample/bivariate/polychoric/processor.rb', line 78

def fd_loglike_cell_b(i, j, m)
  if m==j
     eq12b(a(i),b(m))-eq12b(a(i-1),b(m))
  elsif m==(j-1)
    -eq12b(a(i),b(m))+eq12b(a(i-1),b(m))
  else
    0
  end
=begin          
  if m==j
    Distribution::NormalBivariate.pd_cdf_x(a(i),b(m), rho) - Distribution::NormalBivariate.pd_cdf_x(a(i-1),b(m),rho)
  elsif m==(j-1)
    -Distribution::NormalBivariate.pd_cdf_x(a(i),b(m),rho) + Distribution::NormalBivariate.pd_cdf_x(a(i-1),b(m),rho)
  else
    0
  end
=end          
  
  
end

#fd_loglike_cell_rho(i, j) ⇒ Object

Equation(8) from Olsson(1979)



55
56
57
 
# File 'lib/statsample/bivariate/polychoric/processor.rb', line 55

def fd_loglike_cell_rho(i, j)
  bipdf(i,j) - bipdf(i-1,j) - bipdf(i, j-1) + bipdf(i-1, j-1)
end

#fd_loglike_rho ⇒ Object

First derivate for rho Uses equation (9) from Olsson(1979)



129
130
131
132
133
134
135
136
137
138
139
140
141
142
 
# File 'lib/statsample/bivariate/polychoric/processor.rb', line 129

def fd_loglike_rho
  rho=@rho
  if rho.abs>0.9999
    rho= (rho>0) ? 0.9999 : -0.9999
  end
  total=0
  @nr.times do |i|
    @nc.times do |j|
      pi=pd[i][j] + EPSILON
      total+= (@matrix[i,j].quo(pi))  * (bipdf(i,j)-bipdf(i-1,j)-bipdf(i,j-1)+bipdf(i-1,j-1))  
    end
  end
  total
end

#im_function(t, i, j) ⇒ Object

Returns the derivative correct according to order



229
230
231
232
233
234
235
236
237
238
239
 
# File 'lib/statsample/bivariate/polychoric/processor.rb', line 229

def im_function(t,i,j)
  if t==0
    fd_loglike_cell_rho(i,j)
  elsif t>=1 and t<=@alpha.size
    fd_loglike_cell_a(i,j,t-1)
  elsif t>=@alpha.size+1 and t<=(@alpha.size+@beta.size)
    fd_loglike_cell_b(i,j,t-@alpha.size-1)
  else
    raise "incorrect #{t}"
  end
end

#information_matrix ⇒ Object 



240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
 
# File 'lib/statsample/bivariate/polychoric/processor.rb', line 240

def information_matrix
  total_n=@matrix.total_sum
  vars=@alpha.size+@beta.size+1
  matrix=vars.times.map { vars.times.map {0}}
  vars.times do |m|
    vars.times do |n|
      total=0
      (@nr-1).times do |i|
        (@nc-1).times do |j|
          total+=(1.quo(pd[i][j]+EPSILON)) * im_function(m,i,j) * im_function(n,i,j)
        end
      end
      matrix[m][n]=total_n*total
    end
  end
  m=::Matrix.rows(matrix)
   
end

#loglike ⇒ Object 



22
23
24
25
26
27
28
29
30
31
32
33
34
35
 
# File 'lib/statsample/bivariate/polychoric/processor.rb', line 22

def loglike
  rho=@rho
  if rho.abs>0.9999
    rho= (rho>0) ? 0.9999 : -0.9999
  end
  loglike=0
  @nr.times do |i|
    @nc.times do |j|
      res=pd[i][j]+EPSILON
      loglike+= @matrix[i,j]  * Math::log( res )
    end
  end
  -loglike
end

#pd ⇒ Object

phi_ij for each i and j Uses equation(4) from Olsson(1979)



101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
 
# File 'lib/statsample/bivariate/polychoric/processor.rb', line 101

def pd
  if @pd.nil?
    @pd=@nr.times.collect{ [0] * @nc}
    pc=@nr.times.collect{ [0] * @nc}
    @nr.times do |i|
    @nc.times do |j|
     
      if i==@nr-1 and j==@nc-1
        @pd[i][j]=1.0
      else
        a=(i==@nr-1) ? 100: alpha[i]
        b=(j==@nc-1) ? 100: beta[j]
        #puts "a:#{a} b:#{b}"
        @pd[i][j]=Distribution::BivariateNormal.cdf(a, b, rho)
      end
      pc[i][j] = @pd[i][j]
      @pd[i][j] = @pd[i][j] - pc[i-1][j] if i>0
      @pd[i][j] = @pd[i][j] - pc[i][j-1] if j>0
      @pd[i][j] = @pd[i][j] + pc[i-1][j-1] if (i>0 and j>0)
    end
    end
  end
  @pd
end