Parent

Statsample::Bivariate::Polychoric::Processor

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

Constants

EPSILON

Attributes

alpha[R]
beta[R]
rho[R]
matrix[R]

Public Class Methods

new(alpha,beta,rho,matrix=nil) click to toggle source
    # File lib/statsample/bivariate/polychoric/processor.rb, line 8
 8:         def initialize(alpha,beta,rho,matrix=nil)
 9:           @alpha=alpha
10:           @beta=beta
11:           @matrix=matrix
12:           @nr=@alpha.size+1
13:           @nc=@beta.size+1
14:           @rho=rho
15:           @pd=nil
16:         end

Public Instance Methods

a(i) click to toggle source
    # File lib/statsample/bivariate/polychoric/processor.rb, line 37
37:         def a(i)
38:           raise "Index #{i} should be <= #{@nr-1}" if i>@nr-1
39:           i < 0 ? 100 : (i==@nr-1 ? 100 : alpha[i])
40:         end
b(j) click to toggle source
    # File lib/statsample/bivariate/polychoric/processor.rb, line 41
41:         def b(j)
42:           raise "Index #{j} should be <= #{@nc-1}" if j>@nc-1
43:           j < 0 ? 100 : (j==@nc-1 ? 100 : beta[j])
44:         end
bipdf(i,j) click to toggle source
    # File lib/statsample/bivariate/polychoric/processor.rb, line 18
18:         def bipdf(i,j)
19:            Distribution::NormalBivariate.pdf(a(i), b(j), rho)
20:         end
eq12(u,v) click to toggle source
    # File lib/statsample/bivariate/polychoric/processor.rb, line 46
46:         def eq12(u,v)
47:           Distribution::Normal.pdf(u)*Distribution::Normal.cdf((v-rho*u).quo( Math::sqrt(1-rho**2)))
48:         end
eq12b(u,v) click to toggle source
    # File lib/statsample/bivariate/polychoric/processor.rb, line 50
50:         def eq12b(u,v)
51:           Distribution::Normal.pdf(v) * Distribution::Normal.cdf((u-rho*v).quo( Math::sqrt(1-rho**2)))
52:           
53:         end
fd_loglike_a(k) click to toggle source

First derivative for alpha_k Uses equation (6)

     # File lib/statsample/bivariate/polychoric/processor.rb, line 146
146:         def fd_loglike_a(k)
147:           fd_loglike_a_eq6(k)
148:         end
fd_loglike_a_eq13(k) click to toggle source

Uses equation(13) from Olsson(1979)

     # File lib/statsample/bivariate/polychoric/processor.rb, line 169
169:         def fd_loglike_a_eq13(k)
170:           rho=@rho
171:           if rho.abs>0.9999
172:             rho= (rho>0) ? 0.9999 : 0.9999
173:           end
174:           total=0
175:           a_k=a(k)
176:           @nc.times do |j|
177:             #puts "j: #{j}"
178:             #puts "b #{j} : #{b.call(j)}"
179:             #puts "b #{j-1} : #{b.call(j-1)}"
180:           
181:             e_1=@matrix[k,j].quo(pd[k][j]+EPSILON) - @matrix[k+1,j].quo(pd[k+1][j]+EPSILON)
182:             e_2=Distribution::Normal.pdf(a_k)
183:             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)))
184:             #puts "val #{j}: #{e_1} | #{e_2} | #{e_3}"
185:             total+= e_1*e_2*e_3
186:           end
187:           total
188:         end
fd_loglike_a_eq6(k) click to toggle source

Uses equation (6) from Olsson(1979)

     # File lib/statsample/bivariate/polychoric/processor.rb, line 153
153:         def fd_loglike_a_eq6(k)
154:           rho=@rho
155:           if rho.abs>0.9999
156:             rho= (rho>0) ? 0.9999 : 0.9999
157:           end
158:           total=0
159:           @nr.times do |i|
160:             @nc.times  do |j|
161:               total+=@matrix[i,j].quo(pd[i][j]+EPSILON) * fd_loglike_cell_a(i,j,k)
162:             end
163:           end
164:           total
165:         end
fd_loglike_b(m) click to toggle source

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

     # File lib/statsample/bivariate/polychoric/processor.rb, line 206
206:         def fd_loglike_b(m)
207:           fd_loglike_b_eq14(m)
208:         end
fd_loglike_b_eq14(m) click to toggle source

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

     # File lib/statsample/bivariate/polychoric/processor.rb, line 211
211:         def fd_loglike_b_eq14(m)
212:           rho=@rho
213:           if rho.abs>0.9999
214:             rho= (rho>0) ? 0.9999 : 0.9999
215:           end
216:           total=0
217:           b_m=b(m)
218:           @nr.times do |i|
219:             e_1=@matrix[i,m].quo(pd[i][m]+EPSILON) - @matrix[i,m+1].quo(pd[i][m+1]+EPSILON)
220:             e_2=Distribution::Normal.pdf(b_m)
221:             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)))
222:             #puts "val #{j}: #{e_1} | #{e_2} | #{e_3}"
223:             
224:             total+= e_1*e_2*e_3
225:           end
226:           total
227:         end
fd_loglike_b_eq6(m) click to toggle source

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

     # File lib/statsample/bivariate/polychoric/processor.rb, line 191
191:         def fd_loglike_b_eq6(m)
192:           rho=@rho
193:           if rho.abs>0.9999
194:             rho= (rho>0) ? 0.9999 : 0.9999
195:           end
196:           total=0
197:           @nr.times do |i|
198:             @nc.times  do |j|
199:               total+=@matrix[i,j].quo(pd[i][j]+EPSILON) * fd_loglike_cell_b(i,j,m)
200:             end
201:           end
202:           total
203:         end
fd_loglike_cell_a(i, j, k) click to toggle source

Equation(10) from Olsson(1979)

    # File lib/statsample/bivariate/polychoric/processor.rb, line 59
59:         def fd_loglike_cell_a(i, j, k)
60:           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          
61:           if k==i
62:             eq12(a(k),b(j))-eq12(a(k), b(j-1))
63:           elsif k==(i-1)
64:             -eq12(a(k),b(j))+eq12(a(k), b(j-1))
65:           else
66:             0
67:           end          
68:         end
fd_loglike_cell_b(i, j, m) click to toggle source
    # File lib/statsample/bivariate/polychoric/processor.rb, line 78
78:         def fd_loglike_cell_b(i, j, m)
79:           if m==j
80:              eq12b(a(i),b(m))-eq12b(a(i-1),b(m))
81:           elsif m==(j-1)
82:             -eq12b(a(i),b(m))+eq12b(a(i-1),b(m))
83:           else
84:             0
85:           end
86:                     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          
87:           
88:           
89:         end
fd_loglike_cell_rho(i, j) click to toggle source

Equation(8) from Olsson(1979)

    # File lib/statsample/bivariate/polychoric/processor.rb, line 55
55:         def fd_loglike_cell_rho(i, j)
56:           bipdf(i,j) - bipdf(i-1,j) - bipdf(i, j-1) + bipdf(i-1, j-1)
57:         end
fd_loglike_rho() click to toggle source

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

     # File lib/statsample/bivariate/polychoric/processor.rb, line 129
129:         def fd_loglike_rho
130:           rho=@rho
131:           if rho.abs>0.9999
132:             rho= (rho>0) ? 0.9999 : 0.9999
133:           end
134:           total=0
135:           @nr.times do |i|
136:             @nc.times do |j|
137:               pi=pd[i][j] + EPSILON
138:               total+= (@matrix[i,j].quo(pi))  * (bipdf(i,j)-bipdf(i-1,j)-bipdf(i,j-1)+bipdf(i-1,j-1))  
139:             end
140:           end
141:           total
142:         end
im_function(t,i,j) click to toggle source

Returns the derivative correct according to order

     # File lib/statsample/bivariate/polychoric/processor.rb, line 229
229:         def im_function(t,i,j)
230:           if t==0
231:             fd_loglike_cell_rho(i,j)
232:           elsif t>=1 and t<=@alpha.size
233:             fd_loglike_cell_a(i,j,t-1)
234:           elsif t>=@alpha.size+1 and t<=(@alpha.size+@beta.size)
235:             fd_loglike_cell_b(i,j,t-@alpha.size-1)
236:           else
237:             raise "incorrect #{t}"
238:           end
239:         end
information_matrix() click to toggle source
     # File lib/statsample/bivariate/polychoric/processor.rb, line 240
240:         def information_matrix
241:           total_n=@matrix.total_sum
242:           vars=@alpha.size+@beta.size+1
243:           matrix=vars.times.map { vars.times.map {0}}
244:           vars.times do |m|
245:             vars.times do |n|
246:               total=0
247:               (@nr-1).times do |i|
248:                 (@nc-1).times do |j|
249:                   total+=(1.quo(pd[i][j]+EPSILON)) * im_function(m,i,j) * im_function(n,i,j)
250:                 end
251:               end
252:               matrix[m][n]=total_n*total
253:             end
254:           end
255:           m=::Matrix.rows(matrix)
256:            
257:         end
loglike() click to toggle source
    # File lib/statsample/bivariate/polychoric/processor.rb, line 22
22:         def loglike
23:           rho=@rho
24:           if rho.abs>0.9999
25:             rho= (rho>0) ? 0.9999 : 0.9999
26:           end
27:           loglike=0
28:           @nr.times do |i|
29:             @nc.times do |j|
30:               res=pd[i][j]+EPSILON
31:               loglike+= @matrix[i,j]  * Math::log( res )
32:             end
33:           end
34:           -loglike
35:         end
pd() click to toggle source

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

     # File lib/statsample/bivariate/polychoric/processor.rb, line 101
101:         def pd
102:           if @pd.nil?
103:             @pd=@nr.times.collect{ [0] * @nc}
104:             pc=@nr.times.collect{ [0] * @nc}
105:             @nr.times do |i|
106:             @nc.times do |j|
107:              
108:               if i==@nr-1 and j==@nc-1
109:                 @pd[i][j]=1.0
110:               else
111:                 a=(i==@nr-1) ? 100: alpha[i]
112:                 b=(j==@nc-1) ? 100: beta[j]
113:                 #puts "a:#{a} b:#{b}"
114:                 @pd[i][j]=Distribution::NormalBivariate.cdf(a, b, rho)
115:               end
116:               pc[i][j] = @pd[i][j]
117:               @pd[i][j] = @pd[i][j] - pc[i-1][j] if i>0
118:               @pd[i][j] = @pd[i][j] - pc[i][j-1] if j>0
119:               @pd[i][j] = @pd[i][j] + pc[i-1][j-1] if (i>0 and j>0)
120:             end
121:             end
122:           end
123:           @pd
124:         end

Disabled; run with --debug to generate this.

[Validate]

Generated with the Darkfish Rdoc Generator 1.1.6.