CROSSTABS: Calculate ASE for asymmetric lambda (correctly).

[pspp] / src / language / stats / crosstabs.q

diff --git a/src/language/stats/crosstabs.q b/src/language/stats/crosstabs.q
index 5095b7471551bd57fad15da0db9eb16ff6a8e182..be22e62600c0bb72c21fbba98204196f00e6963d 100644 (file)
--- a/src/language/stats/crosstabs.q
+++ b/src/language/stats/crosstabs.q
@@ -16,12 +16,10 @@
 
 /* FIXME:
 
-   - T values for Spearman's R and Pearson's R are wrong.
    - How to calculate significance of symmetric and directional measures?
-   - Asymmetric ASEs and T values for lambda are wrong.
-   - ASE of Goodman and Kruskal's tau is not calculated.
-   - ASE of symmetric somers' d is wrong.
-   - Approx. T of uncertainty coefficient is wrong.
+   - How to calculate ASE for symmetric Somers ' d?
+   - How to calculate ASE for Goodman and Kruskal's tau?
+   - How to calculate approx. T of symmetric uncertainty coefficient?
 
 */
 
@@ -2265,12 +2263,12 @@ calc_chisq (struct pivot_table *pt,
     }
 }
 
-/* Calculate the value of Pearson's r.  r is stored into R, ase_1 into
-   ASE_1, and ase_0 into ASE_0.  The row and column values must be
+/* Calculate the value of Pearson's r.  r is stored into R, its T value into
+   T, and standard error into ERROR.  The row and column values must be
    passed in PT and Y. */
 static void
 calc_r (struct pivot_table *pt,
-        double *PT, double *Y, double *r, double *ase_0, double *ase_1)
+        double *PT, double *Y, double *r, double *t, double *error)
 {
   double SX, SY, S, T;
   double Xbar, Ybar;
@@ -2308,7 +2306,7 @@ calc_r (struct pivot_table *pt,
   SY = sum_Y2c - pow2 (sum_Yc) / pt->total;
   T = sqrt (SX * SY);
   *r = S / T;
-  *ase_0 = sqrt ((sum_X2Y2f - pow2 (sum_XYf) / pt->total) / (sum_X2r * sum_Y2c));
+  *t = *r / sqrt (1 - pow2 (*r)) * sqrt (pt->total - 2);
 
   {
     double s, c, y, t;
@@ -2329,7 +2327,7 @@ calc_r (struct pivot_table *pt,
          c = (t - s) - y;
          s = t;
        }
-    *ase_1 = sqrt (s) / (T * T);
+    *error = sqrt (s) / (T * T);
   }
 }
 
@@ -2536,7 +2534,7 @@ calc_symmetric (struct crosstabs_proc *proc, struct pivot_table *pt,
       if (proc->statistics & (1u << CRS_ST_D))
        {
          somers_d_v[0] = (P - Q) / (.5 * (Dc + Dr));
-         somers_d_ase[0] = 2. * btau_var / (Dr + Dc) * sqrt (Dr * Dc);
+         somers_d_ase[0] = SYSMIS;
          somers_d_t[0] = (somers_d_v[0]
                           / (4 / (Dc + Dr)
                              * sqrt (ctau_cum - pow2 (P - Q) / pt->total)));
@@ -2596,13 +2594,11 @@ calc_symmetric (struct crosstabs_proc *proc, struct pivot_table *pt,
       }
 
       calc_r (pt, R, C, &v[6], &t[6], &ase[6]);
-      t[6] = v[6] / t[6];
 
       free (R);
       free (C);
 
       calc_r (pt, (double *) pt->rows, (double *) pt->cols, &v[7], &t[7], &ase[7]);
-      t[7] = v[7] / t[7];
     }
 
   /* Cohen's kappa. */
@@ -2805,19 +2801,14 @@ calc_directional (struct crosstabs_proc *proc, struct pivot_table *pt,
 
       /* ASE1 for Y given PT. */
       {
-       double accum;
+        double accum;
 
-       for (accum = 0., i = 0; i < pt->n_rows; i++)
-         for (j = 0; j < pt->n_cols; j++)
-           {
-             const int deltaj = j == cm_index;
-             accum += (pt->mat[j + i * pt->n_cols]
-                       * pow2 ((j == fim_index[i])
-                              - deltaj
-                              + v[0] * deltaj));
-           }
-
-       ase[2] = sqrt (accum - pt->total * v[0]) / (pt->total - cm);
+        accum = 0.;
+       for (i = 0; i < pt->n_rows; i++)
+          if (cm_index == fim_index[i])
+            accum += fim[i];
+        ase[2] = sqrt ((pt->total - sum_fim) * (sum_fim + cm - 2. * accum)
+                       / pow3 (pt->total - cm));
       }
 
       /* ASE0 for Y given PT. */
@@ -2833,19 +2824,14 @@ calc_directional (struct crosstabs_proc *proc, struct pivot_table *pt,
 
       /* ASE1 for PT given Y. */
       {
-       double accum;
+        double accum;
 
-       for (accum = 0., i = 0; i < pt->n_rows; i++)
-         for (j = 0; j < pt->n_cols; j++)
-           {
-             const int deltaj = i == rm_index;
-             accum += (pt->mat[j + i * pt->n_cols]
-                       * pow2 ((i == fmj_index[j])
-                              - deltaj
-                              + v[0] * deltaj));
-           }
-
-       ase[1] = sqrt (accum - pt->total * v[0]) / (pt->total - rm);
+        accum = 0.;
+       for (j = 0; j < pt->n_cols; j++)
+          if (rm_index == fmj_index[j])
+            accum += fmj[j];
+        ase[1] = sqrt ((pt->total - sum_fmj) * (sum_fmj + rm - 2. * accum)
+                       / pow3 (pt->total - rm));
       }
 
       /* ASE0 for PT given Y. */
@@ -2874,7 +2860,7 @@ calc_directional (struct crosstabs_proc *proc, struct pivot_table *pt,
                         * pow2 (temp0 + (v[0] - 1.) * temp1));
            }
        ase[0] = sqrt (accum1 - 4. * pt->total * v[0] * v[0]) / (2. * pt->total - rm - cm);
-       t[0] = v[0] / (sqrt (accum0 - pow2 ((sum_fim + sum_fmj - cm - rm) / pt->total))
+       t[0] = v[0] / (sqrt (accum0 - pow2 (sum_fim + sum_fmj - cm - rm) / pt->total)
                       / (2. * pt->total - rm - cm));
       }
 
@@ -2951,8 +2937,7 @@ calc_directional (struct crosstabs_proc *proc, struct pivot_table *pt,
 
       v[5] = 2. * ((UX + UY - UXY) / (UX + UY));
       ase[5] = (2. / (pt->total * pow2 (UX + UY))) * sqrt (ase1_sym);
-      t[5] = v[5] / ((2. / (pt->total * (UX + UY)))
-                    * sqrt (P - pow2 (UX + UY - UXY) / pt->total));
+      t[5] = SYSMIS;
 
       v[6] = (UX + UY - UXY) / UX;
       ase[6] = sqrt (ase1_xy) / (pt->total * UX * UX);