polychoric.rb in statsample-0.10.0

- old
+ new

@@ -8,22 +8,33 @@
     end
     
     # Polychoric correlation matrix.
     # Order of rows and columns depends on Dataset#fields order
     def self.polychoric_correlation_matrix(ds)
-      ds.collect_matrix do |row,col|
+      cache={}
+      matrix=ds.collect_matrix do |row,col|
         if row==col
           1.0
         else
           begin
-            polychoric(ds[row],ds[col])
+            if cache[[col,row]].nil?
+              poly=polychoric(ds[row],ds[col])
+              cache[[row,col]]=poly
+              poly
+            else
+              cache[[col,row]]
+            end
           rescue RuntimeError
             nil
           end
         end
       end
+      matrix.extend CovariateMatrix
+      matrix.fields=ds.fields
+      matrix
     end
+    
     # = Polychoric correlation.
     #
     # The <em>polychoric</em> correlation is a measure of 
     # bivariate association arising when both observed variates
     # are ordered, categorical variables that result from polychotomizing
@@ -81,10 +92,11 @@
       # See http://rb-gsl.rubyforge.org/min.html for reference.  
       attr_accessor :minimizer_type_joint
       
       
       # Method of calculation of polychoric series. 
+      # <tt>:two_step</tt> used by default.
       # 
       # :two_step:: two-step ML, based on code by Gegenfurtner(1992).
       # :polychoric_series:: polychoric series estimate, using 
       #                      algorithm AS87 by Martinson and Hamdan (1975).
       # :joint::             one-step ML, based on R package 'polycor'
@@ -105,11 +117,11 @@
       METHOD=:two_step
       MAX_ITERATIONS=300
       EPSILON=1e-6
       MINIMIZER_TYPE_TWO_STEP="brent"
       MINIMIZER_TYPE_JOINT="nmsimplex"
-      def new_with_vectors(v1,v2)
+      def self.new_with_vectors(v1,v2)
         Polychoric.new(Crosstab.new(v1,v2).to_matrix)
       end
       # Params:
       # * matrix: Contingence table
       # * opts: Any attribute
@@ -247,18 +259,19 @@
             else
               a=(i==@nr-1) ? 100: alpha[i]
               b=(j==@nc-1) ? 100: beta[j]
               #puts "a:#{a} b:#{b}"
               pd[i][j]=Distribution::NormalBivariate.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)
             res= pd[i][j]
              #puts "i:#{i} | j:#{j} | ac: #{sprintf("%0.4f", pc[i][j])} | pd: #{sprintf("%0.4f", pd[i][j])} | res:#{sprintf("%0.4f", res)}"
-          if (res==0)
+          if (res<=0)
            #    puts "Correccion"
             res=1e-16
           end
             loglike+= @matrix[i,j]  * Math::log( res )
           }
@@ -326,10 +339,10 @@
         @log="Minimizing using GSL Brent method\n"
         min=Minimization::Brent.new(-0.9999,0.9999,f)
         min.epsilon=@epsilon
         min.expected=0
         min.iterate
-        @log+=min.log
+        @log+=min.log.to_table.to_s
         @r=min.x_minimum
         @loglike_model=-min.f_minimum
         puts @log if @debug
         
       end