pintos-os.org Git - pspp/blob - src/language/stats/factor.c

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 2009, 2010, 2011 Free Software Foundation, Inc.
   3
   4    This program is free software: you can redistribute it and/or modify
   5    it under the terms of the GNU General Public License as published by
   6    the Free Software Foundation, either version 3 of the License, or
   7    (at your option) any later version.
   8
   9    This program is distributed in the hope that it will be useful,
  10    but WITHOUT ANY WARRANTY; without even the implied warranty of
  11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12    GNU General Public License for more details.
  13
  14    You should have received a copy of the GNU General Public License
  15    along with this program.  If not, see <http://www.gnu.org/licenses/>. */
  16
  17 #include <config.h>
  18
  19 #include <gsl/gsl_vector.h>
  20 #include <gsl/gsl_linalg.h>
  21 #include <gsl/gsl_matrix.h>
  22 #include <gsl/gsl_eigen.h>
  23 #include <gsl/gsl_blas.h>
  24 #include <gsl/gsl_sort_vector.h>
  25 #include <gsl/gsl_cdf.h>
  26
  27 #include "data/casegrouper.h"
  28 #include "data/casereader.h"
  29 #include "data/casewriter.h"
  30 #include "data/dataset.h"
  31 #include "data/dictionary.h"
  32 #include "data/format.h"
  33 #include "data/subcase.h"
  34 #include "language/command.h"
  35 #include "language/lexer/lexer.h"
  36 #include "language/lexer/value-parser.h"
  37 #include "language/lexer/variable-parser.h"
  38 #include "libpspp/message.h"
  39 #include "libpspp/misc.h"
  40 #include "math/correlation.h"
  41 #include "math/covariance.h"
  42 #include "math/moments.h"
  43 #include "output/chart-item.h"
  44 #include "output/charts/scree.h"
  45 #include "output/tab.h"
  46
  47 #include "gettext.h"
  48 #define _(msgid) gettext (msgid)
  49 #define N_(msgid) msgid
  50
  51 enum method
  52   {
  53     METHOD_CORR,
  54     METHOD_COV
  55   };
  56
  57 enum missing_type
  58   {
  59     MISS_LISTWISE,
  60     MISS_PAIRWISE,
  61     MISS_MEANSUB,
  62   };
  63
  64 enum extraction_method
  65   {
  66     EXTRACTION_PC,
  67     EXTRACTION_PAF,
  68   };
  69
  70 enum plot_opts
  71   {
  72     PLOT_SCREE = 0x0001,
  73     PLOT_ROTATION = 0x0002
  74   };
  75
  76 enum print_opts
  77   {
  78     PRINT_UNIVARIATE  = 0x0001,
  79     PRINT_DETERMINANT = 0x0002,
  80     PRINT_INV         = 0x0004,
  81     PRINT_AIC         = 0x0008,
  82     PRINT_SIG         = 0x0010,
  83     PRINT_COVARIANCE  = 0x0020,
  84     PRINT_CORRELATION = 0x0040,
  85     PRINT_ROTATION    = 0x0080,
  86     PRINT_EXTRACTION  = 0x0100,
  87     PRINT_INITIAL     = 0x0200,
  88     PRINT_KMO         = 0x0400,
  89     PRINT_REPR        = 0x0800,
  90     PRINT_FSCORE      = 0x1000
  91   };
  92
  93 enum rotation_type
  94   {
  95     ROT_VARIMAX = 0,
  96     ROT_EQUAMAX,
  97     ROT_QUARTIMAX,
  98     ROT_NONE
  99   };
 100
 101 typedef void (*rotation_coefficients) (double *x, double *y,
 102                                     double a, double b, double c, double d,
 103                                     const gsl_matrix *loadings );
 104
 105
 106 static void
 107 varimax_coefficients (double *x, double *y,
 108                       double a, double b, double c, double d,
 109                       const gsl_matrix *loadings )
 110 {
 111   *x = d - 2 * a * b / loadings->size1;
 112   *y = c - (a * a - b * b) / loadings->size1;
 113 }
 114
 115 static void
 116 equamax_coefficients (double *x, double *y,
 117                       double a, double b, double c, double d,
 118                       const gsl_matrix *loadings )
 119 {
 120   *x = d - loadings->size2 * a * b / loadings->size1;
 121   *y = c - loadings->size2 * (a * a - b * b) / (2 * loadings->size1);
 122 }
 123
 124 static void
 125 quartimax_coefficients (double *x, double *y,
 126                       double a UNUSED, double b UNUSED, double c, double d,
 127                       const gsl_matrix *loadings UNUSED)
 128 {
 129   *x = d ;
 130   *y = c ;
 131 }
 132
 133 static const rotation_coefficients rotation_coeff[3] = {
 134   varimax_coefficients,
 135   equamax_coefficients,
 136   quartimax_coefficients
 137 };
 138
 139
 140 struct cmd_factor
 141 {
 142   size_t n_vars;
 143   const struct variable **vars;
 144
 145   const struct variable *wv;
 146
 147   enum method method;
 148   enum missing_type missing_type;
 149   enum mv_class exclude;
 150   enum print_opts print;
 151   enum extraction_method extraction;
 152   enum plot_opts plot;
 153   enum rotation_type rotation;
 154
 155   /* Extraction Criteria */
 156   int n_factors;
 157   double min_eigen;
 158   double econverge;
 159   int iterations;
 160
 161   double rconverge;
 162
 163   /* Format */
 164   double blank;
 165   bool sort;
 166 };
 167
 168 struct idata
 169 {
 170   /* Intermediate values used in calculation */
 171
 172   const gsl_matrix *corr ;  /* The correlation matrix */
 173   gsl_matrix *cov ;         /* The covariance matrix */
 174   const gsl_matrix *n ;     /* Matrix of number of samples */
 175
 176   gsl_vector *eval ;  /* The eigenvalues */
 177   gsl_matrix *evec ;  /* The eigenvectors */
 178
 179   int n_extractions;
 180
 181   gsl_vector *msr ;  /* Multiple Squared Regressions */
 182
 183   double detR;  /* The determinant of the correlation matrix */
 184 };
 185
 186 static struct idata *
 187 idata_alloc (size_t n_vars)
 188 {
 189   struct idata *id = xzalloc (sizeof (*id));
 190
 191   id->n_extractions = 0;
 192   id->msr = gsl_vector_alloc (n_vars);
 193
 194   id->eval = gsl_vector_alloc (n_vars);
 195   id->evec = gsl_matrix_alloc (n_vars, n_vars);
 196
 197   return id;
 198 }
 199
 200 static void
 201 idata_free (struct idata *id)
 202 {
 203   gsl_vector_free (id->msr);
 204   gsl_vector_free (id->eval);
 205   gsl_matrix_free (id->evec);
 206   if (id->cov != NULL)
 207     gsl_matrix_free (id->cov);
 208
 209   free (id);
 210 }
 211
 212
 213 static gsl_matrix *
 214 anti_image (const gsl_matrix *m)
 215 {
 216   int i, j;
 217   gsl_matrix *a;
 218   assert (m->size1 == m->size2);
 219
 220   a = gsl_matrix_alloc (m->size1, m->size2);
 221
 222   for (i = 0; i < m->size1; ++i)
 223     {
 224       for (j = 0; j < m->size2; ++j)
 225         {
 226           double *p = gsl_matrix_ptr (a, i, j);
 227           *p = gsl_matrix_get (m, i, j);
 228           *p /= gsl_matrix_get (m, i, i);
 229           *p /= gsl_matrix_get (m, j, j);
 230         }
 231     }
 232
 233   return a;
 234 }
 235
 236
 237 /* Return the sum of all the elements excluding row N */
 238 static double
 239 ssq_od_n (const gsl_matrix *m, int n)
 240 {
 241   int i, j;
 242   double ss = 0;
 243   assert (m->size1 == m->size2);
 244
 245   assert (n < m->size1);
 246
 247   for (i = 0; i < m->size1; ++i)
 248     {
 249       if (i == n ) continue;
 250       for (j = 0; j < m->size2; ++j)
 251         {
 252           ss += pow2 (gsl_matrix_get (m, i, j));
 253         }
 254     }
 255
 256   return ss;
 257 }
 258
 259
 260
 261 #if 0
 262 static void
 263 dump_matrix (const gsl_matrix *m)
 264 {
 265   size_t i, j;
 266
 267   for (i = 0 ; i < m->size1; ++i)
 268     {
 269       for (j = 0 ; j < m->size2; ++j)
 270         printf ("%02f ", gsl_matrix_get (m, i, j));
 271       printf ("\n");
 272     }
 273 }
 274
 275 static void
 276 dump_matrix_permute (const gsl_matrix *m, const gsl_permutation *p)
 277 {
 278   size_t i, j;
 279
 280   for (i = 0 ; i < m->size1; ++i)
 281     {
 282       for (j = 0 ; j < m->size2; ++j)
 283         printf ("%02f ", gsl_matrix_get (m, gsl_permutation_get (p, i), j));
 284       printf ("\n");
 285     }
 286 }
 287
 288
 289 static void
 290 dump_vector (const gsl_vector *v)
 291 {
 292   size_t i;
 293   for (i = 0 ; i < v->size; ++i)
 294     {
 295       printf ("%02f\n", gsl_vector_get (v, i));
 296     }
 297   printf ("\n");
 298 }
 299 #endif
 300
 301
 302 static int
 303 n_extracted_factors (const struct cmd_factor *factor, struct idata *idata)
 304 {
 305   int i;
 306
 307   /* If there is a cached value, then return that. */
 308   if ( idata->n_extractions != 0)
 309     return idata->n_extractions;
 310
 311   /* Otherwise, if the number of factors has been explicitly requested,
 312      use that. */
 313   if (factor->n_factors > 0)
 314     {
 315       idata->n_extractions = factor->n_factors;
 316       goto finish;
 317     }
 318
 319   /* Use the MIN_EIGEN setting. */
 320   for (i = 0 ; i < idata->eval->size; ++i)
 321     {
 322       double evali = fabs (gsl_vector_get (idata->eval, i));
 323
 324       idata->n_extractions = i;
 325
 326       if (evali < factor->min_eigen)
 327         goto finish;
 328     }
 329
 330  finish:
 331   return idata->n_extractions;
 332 }
 333
 334
 335 /* Returns a newly allocated matrix identical to M.
 336    It it the callers responsibility to free the returned value.
 337 */
 338 static gsl_matrix *
 339 matrix_dup (const gsl_matrix *m)
 340 {
 341   gsl_matrix *n =  gsl_matrix_alloc (m->size1, m->size2);
 342
 343   gsl_matrix_memcpy (n, m);
 344
 345   return n;
 346 }
 347
 348
 349 struct smr_workspace
 350 {
 351   /* Copy of the subject */
 352   gsl_matrix *m;
 353
 354   gsl_matrix *inverse;
 355
 356   gsl_permutation *perm;
 357
 358   gsl_matrix *result1;
 359   gsl_matrix *result2;
 360 };
 361
 362
 363 static struct smr_workspace *ws_create (const gsl_matrix *input)
 364 {
 365   struct smr_workspace *ws = xmalloc (sizeof (*ws));
 366
 367   ws->m = gsl_matrix_alloc (input->size1, input->size2);
 368   ws->inverse = gsl_matrix_calloc (input->size1 - 1, input->size2 - 1);
 369   ws->perm = gsl_permutation_alloc (input->size1 - 1);
 370   ws->result1 = gsl_matrix_calloc (input->size1 - 1, 1);
 371   ws->result2 = gsl_matrix_calloc (1, 1);
 372
 373   return ws;
 374 }
 375
 376 static void
 377 ws_destroy (struct smr_workspace *ws)
 378 {
 379   gsl_matrix_free (ws->result2);
 380   gsl_matrix_free (ws->result1);
 381   gsl_permutation_free (ws->perm);
 382   gsl_matrix_free (ws->inverse);
 383   gsl_matrix_free (ws->m);
 384
 385   free (ws);
 386 }
 387
 388
 389 /*
 390    Return the square of the regression coefficient for VAR regressed against all other variables.
 391  */
 392 static double
 393 squared_multiple_correlation (const gsl_matrix *corr, int var, struct smr_workspace *ws)
 394 {
 395   /* For an explanation of what this is doing, see
 396      http://www.visualstatistics.net/Visual%20Statistics%20Multimedia/multiple_regression_analysis.htm
 397   */
 398
 399   int signum = 0;
 400   gsl_matrix_view rxx;
 401
 402   gsl_matrix_memcpy (ws->m, corr);
 403
 404   gsl_matrix_swap_rows (ws->m, 0, var);
 405   gsl_matrix_swap_columns (ws->m, 0, var);
 406
 407   rxx = gsl_matrix_submatrix (ws->m, 1, 1, ws->m->size1 - 1, ws->m->size1 - 1);
 408
 409   gsl_linalg_LU_decomp (&rxx.matrix, ws->perm, &signum);
 410
 411   gsl_linalg_LU_invert (&rxx.matrix, ws->perm, ws->inverse);
 412
 413   {
 414     gsl_matrix_const_view rxy = gsl_matrix_const_submatrix (ws->m, 1, 0, ws->m->size1 - 1, 1);
 415     gsl_matrix_const_view ryx = gsl_matrix_const_submatrix (ws->m, 0, 1, 1, ws->m->size1 - 1);
 416
 417     gsl_blas_dgemm (CblasNoTrans,  CblasNoTrans,
 418                     1.0, ws->inverse, &rxy.matrix, 0.0, ws->result1);
 419
 420     gsl_blas_dgemm (CblasNoTrans,  CblasNoTrans,
 421                     1.0, &ryx.matrix, ws->result1, 0.0, ws->result2);
 422   }
 423
 424   return gsl_matrix_get (ws->result2, 0, 0);
 425 }
 426
 427
 428
 429 static double the_communality (const gsl_matrix *evec, const gsl_vector *eval, int n, int n_factors);
 430
 431
 432 struct factor_matrix_workspace
 433 {
 434   size_t n_factors;
 435   gsl_eigen_symmv_workspace *eigen_ws;
 436
 437   gsl_vector *eval ;
 438   gsl_matrix *evec ;
 439
 440   gsl_matrix *gamma ;
 441
 442   gsl_matrix *r;
 443 };
 444
 445 static struct factor_matrix_workspace *
 446 factor_matrix_workspace_alloc (size_t n, size_t nf)
 447 {
 448   struct factor_matrix_workspace *ws = xmalloc (sizeof (*ws));
 449
 450   ws->n_factors = nf;
 451   ws->gamma = gsl_matrix_calloc (nf, nf);
 452   ws->eigen_ws = gsl_eigen_symmv_alloc (n);
 453   ws->eval = gsl_vector_alloc (n);
 454   ws->evec = gsl_matrix_alloc (n, n);
 455   ws->r  = gsl_matrix_alloc (n, n);
 456
 457   return ws;
 458 }
 459
 460 static void
 461 factor_matrix_workspace_free (struct factor_matrix_workspace *ws)
 462 {
 463   gsl_eigen_symmv_free (ws->eigen_ws);
 464   gsl_vector_free (ws->eval);
 465   gsl_matrix_free (ws->evec);
 466   gsl_matrix_free (ws->gamma);
 467   gsl_matrix_free (ws->r);
 468   free (ws);
 469 }
 470
 471 /*
 472   Shift P left by OFFSET places, and overwrite TARGET
 473   with the shifted result.
 474   Positions in TARGET less than OFFSET are unchanged.
 475 */
 476 static void
 477 perm_shift_apply (gsl_permutation *target, const gsl_permutation *p,
 478                   size_t offset)
 479 {
 480   size_t i;
 481   assert (target->size == p->size);
 482   assert (offset <= target->size);
 483
 484   for (i = 0; i < target->size - offset; ++i)
 485     {
 486       target->data[i] = p->data [i + offset];
 487     }
 488 }
 489
 490
 491 /*
 492    Indirectly sort the rows of matrix INPUT, storing the sort order in PERM.
 493    The sort criteria are as follows:
 494
 495    Rows are sorted on the first column, until the absolute value of an
 496    element in a subsequent column  is greater than that of the first
 497    column.  Thereafter, rows will be sorted on the second column,
 498    until the absolute value of an element in a subsequent column
 499    exceeds that of the second column ...
 500 */
 501 static void
 502 sort_matrix_indirect (const gsl_matrix *input, gsl_permutation *perm)
 503 {
 504   const size_t n = perm->size;
 505   const size_t m = input->size2;
 506   int i, j;
 507   gsl_matrix *mat ;
 508   int column_n = 0;
 509   int row_n = 0;
 510   gsl_permutation *p;
 511
 512   assert (perm->size == input->size1);
 513
 514   p = gsl_permutation_alloc (n);
 515
 516   /* Copy INPUT into MAT, discarding the sign */
 517   mat = gsl_matrix_alloc (n, m);
 518   for (i = 0 ; i < mat->size1; ++i)
 519     {
 520       for (j = 0 ; j < mat->size2; ++j)
 521         {
 522           double x = gsl_matrix_get (input, i, j);
 523           gsl_matrix_set (mat, i, j, fabs (x));
 524         }
 525     }
 526
 527   while (column_n < m && row_n < n)
 528     {
 529       gsl_vector_const_view columni = gsl_matrix_const_column (mat, column_n);
 530       gsl_sort_vector_index (p, &columni.vector);
 531
 532       for (i = 0 ; i < n; ++i)
 533         {
 534           gsl_vector_view row = gsl_matrix_row (mat, p->data[n - 1 - i]);
 535           size_t maxindex = gsl_vector_max_index (&row.vector);
 536
 537           if ( maxindex > column_n )
 538             break;
 539
 540           /* All subsequent elements of this row, are of no interest.
 541              So set them all to a highly negative value */
 542           for (j = column_n + 1; j < row.vector.size ; ++j)
 543             gsl_vector_set (&row.vector, j, -DBL_MAX);
 544         }
 545
 546       perm_shift_apply (perm, p, row_n);
 547       row_n += i;
 548
 549       column_n++;
 550     }
 551
 552   gsl_permutation_free (p);
 553   gsl_matrix_free (mat);
 554
 555   assert ( 0 == gsl_permutation_valid (perm));
 556
 557   /* We want the biggest value to be first */
 558   gsl_permutation_reverse (perm);
 559 }
 560
 561
 562 static void
 563 drot_go (double phi, double *l0, double *l1)
 564 {
 565   double r0 = cos (phi) * *l0 + sin (phi) * *l1;
 566   double r1 = - sin (phi) * *l0 + cos (phi) * *l1;
 567
 568   *l0 = r0;
 569   *l1 = r1;
 570 }
 571
 572
 573 static gsl_matrix *
 574 clone_matrix (const gsl_matrix *m)
 575 {
 576   int j, k;
 577   gsl_matrix *c = gsl_matrix_calloc (m->size1, m->size2);
 578
 579   for (j = 0 ; j < c->size1; ++j)
 580     {
 581       for (k = 0 ; k < c->size2; ++k)
 582         {
 583           const double *v = gsl_matrix_const_ptr (m, j, k);
 584           gsl_matrix_set (c, j, k, *v);
 585         }
 586     }
 587
 588   return c;
 589 }
 590
 591
 592 static double
 593 initial_sv (const gsl_matrix *fm)
 594 {
 595   int j, k;
 596
 597   double sv = 0.0;
 598   for (j = 0 ; j < fm->size2; ++j)
 599     {
 600       double l4s = 0;
 601       double l2s = 0;
 602
 603       for (k = j + 1 ; k < fm->size2; ++k)
 604         {
 605           double lambda = gsl_matrix_get (fm, k, j);
 606           double lambda_sq = lambda * lambda;
 607           double lambda_4 = lambda_sq * lambda_sq;
 608
 609           l4s += lambda_4;
 610           l2s += lambda_sq;
 611         }
 612       sv += ( fm->size1 * l4s - (l2s * l2s) ) / (fm->size1 * fm->size1 );
 613     }
 614   return sv;
 615 }
 616
 617 static void
 618 rotate (const struct cmd_factor *cf, const gsl_matrix *unrot,
 619         const gsl_vector *communalities,
 620         gsl_matrix *result,
 621         gsl_vector *rotated_loadings
 622         )
 623 {
 624   int j, k;
 625   int i;
 626   double prev_sv;
 627
 628   /* First get a normalised version of UNROT */
 629   gsl_matrix *normalised = gsl_matrix_calloc (unrot->size1, unrot->size2);
 630   gsl_matrix *h_sqrt = gsl_matrix_calloc (communalities->size, communalities->size);
 631   gsl_matrix *h_sqrt_inv ;
 632
 633   /* H is the diagonal matrix containing the absolute values of the communalities */
 634   for (i = 0 ; i < communalities->size ; ++i)
 635     {
 636       double *ptr = gsl_matrix_ptr (h_sqrt, i, i);
 637       *ptr = fabs (gsl_vector_get (communalities, i));
 638     }
 639
 640   /* Take the square root of the communalities */
 641   gsl_linalg_cholesky_decomp (h_sqrt);
 642
 643
 644   /* Save a copy of h_sqrt and invert it */
 645   h_sqrt_inv = clone_matrix (h_sqrt);
 646   gsl_linalg_cholesky_decomp (h_sqrt_inv);
 647   gsl_linalg_cholesky_invert (h_sqrt_inv);
 648
 649   /* normalised vertion is H^{1/2} x UNROT */
 650   gsl_blas_dgemm (CblasNoTrans,  CblasNoTrans, 1.0, h_sqrt_inv, unrot, 0.0, normalised);
 651
 652   gsl_matrix_free (h_sqrt_inv);
 653
 654
 655   /* Now perform the rotation iterations */
 656
 657   prev_sv = initial_sv (normalised);
 658   for (i = 0 ; i < cf->iterations ; ++i)
 659     {
 660       double sv = 0.0;
 661       for (j = 0 ; j < normalised->size2; ++j)
 662         {
 663           /* These variables relate to the convergence criterium */
 664           double l4s = 0;
 665           double l2s = 0;
 666
 667           for (k = j + 1 ; k < normalised->size2; ++k)
 668             {
 669               int p;
 670               double a = 0.0;
 671               double b = 0.0;
 672               double c = 0.0;
 673               double d = 0.0;
 674               double x, y;
 675               double phi;
 676
 677               for (p = 0; p < normalised->size1; ++p)
 678                 {
 679                   double jv = gsl_matrix_get (normalised, p, j);
 680                   double kv = gsl_matrix_get (normalised, p, k);
 681
 682                   double u = jv * jv - kv * kv;
 683                   double v = 2 * jv * kv;
 684                   a += u;
 685                   b += v;
 686                   c +=  u * u - v * v;
 687                   d += 2 * u * v;
 688                 }
 689
 690               rotation_coeff [cf->rotation] (&x, &y, a, b, c, d, normalised);
 691
 692               phi = atan2 (x,  y) / 4.0 ;
 693
 694               /* Don't bother rotating if the angle is small */
 695               if ( fabs (sin (phi) ) <= pow (10.0, -15.0))
 696                   continue;
 697
 698               for (p = 0; p < normalised->size1; ++p)
 699                 {
 700                   double *lambda0 = gsl_matrix_ptr (normalised, p, j);
 701                   double *lambda1 = gsl_matrix_ptr (normalised, p, k);
 702                   drot_go (phi, lambda0, lambda1);
 703                 }
 704
 705               /* Calculate the convergence criterium */
 706               {
 707                 double lambda = gsl_matrix_get (normalised, k, j);
 708                 double lambda_sq = lambda * lambda;
 709                 double lambda_4 = lambda_sq * lambda_sq;
 710
 711                 l4s += lambda_4;
 712                 l2s += lambda_sq;
 713               }
 714             }
 715           sv += ( normalised->size1 * l4s - (l2s * l2s) ) / (normalised->size1 * normalised->size1 );
 716         }
 717
 718       if ( fabs (sv - prev_sv) <= cf->rconverge)
 719         break;
 720
 721       prev_sv = sv;
 722     }
 723
 724   gsl_blas_dgemm (CblasNoTrans,  CblasNoTrans, 1.0,
 725                   h_sqrt, normalised,  0.0,   result);
 726
 727   gsl_matrix_free (h_sqrt);
 728
 729
 730   /* reflect negative sums and populate the rotated loadings vector*/
 731   for (i = 0 ; i < result->size2; ++i)
 732     {
 733       double ssq = 0.0;
 734       double sum = 0.0;
 735       for (j = 0 ; j < result->size1; ++j)
 736         {
 737           double s = gsl_matrix_get (result, j, i);
 738           ssq += s * s;
 739           sum += gsl_matrix_get (result, j, i);
 740         }
 741
 742       gsl_vector_set (rotated_loadings, i, ssq);
 743
 744       if ( sum < 0 )
 745         for (j = 0 ; j < result->size1; ++j)
 746           {
 747             double *lambda = gsl_matrix_ptr (result, j, i);
 748             *lambda = - *lambda;
 749           }
 750     }
 751 }
 752
 753
 754 /*
 755   Get an approximation for the factor matrix into FACTORS, and the communalities into COMMUNALITIES.
 756   R is the matrix to be analysed.
 757   WS is a pointer to a structure which must have been initialised with factor_matrix_workspace_init.
 758  */
 759 static void
 760 iterate_factor_matrix (const gsl_matrix *r, gsl_vector *communalities, gsl_matrix *factors,
 761                        struct factor_matrix_workspace *ws)
 762 {
 763   size_t i;
 764   gsl_matrix_view mv ;
 765
 766   assert (r->size1 == r->size2);
 767   assert (r->size1 == communalities->size);
 768
 769   assert (factors->size1 == r->size1);
 770   assert (factors->size2 == ws->n_factors);
 771
 772   gsl_matrix_memcpy (ws->r, r);
 773
 774   /* Apply Communalities to diagonal of correlation matrix */
 775   for (i = 0 ; i < communalities->size ; ++i)
 776     {
 777       double *x = gsl_matrix_ptr (ws->r, i, i);
 778       *x = gsl_vector_get (communalities, i);
 779     }
 780
 781   gsl_eigen_symmv (ws->r, ws->eval, ws->evec, ws->eigen_ws);
 782
 783   mv = gsl_matrix_submatrix (ws->evec, 0, 0, ws->evec->size1, ws->n_factors);
 784
 785   /* Gamma is the diagonal matrix containing the absolute values of the eigenvalues */
 786   for (i = 0 ; i < ws->n_factors ; ++i)
 787     {
 788       double *ptr = gsl_matrix_ptr (ws->gamma, i, i);
 789       *ptr = fabs (gsl_vector_get (ws->eval, i));
 790     }
 791
 792   /* Take the square root of gamma */
 793   gsl_linalg_cholesky_decomp (ws->gamma);
 794
 795   gsl_blas_dgemm (CblasNoTrans,  CblasNoTrans, 1.0, &mv.matrix, ws->gamma, 0.0, factors);
 796
 797   for (i = 0 ; i < r->size1 ; ++i)
 798     {
 799       double h = the_communality (ws->evec, ws->eval, i, ws->n_factors);
 800       gsl_vector_set (communalities, i, h);
 801     }
 802 }
 803
 804
 805
 806 static bool run_factor (struct dataset *ds, const struct cmd_factor *factor);
 807
 808
 809 int
 810 cmd_factor (struct lexer *lexer, struct dataset *ds)
 811 {
 812   bool extraction_seen = false;
 813   const struct dictionary *dict = dataset_dict (ds);
 814
 815   struct cmd_factor factor;
 816   factor.n_vars = 0;
 817   factor.vars = NULL;
 818   factor.method = METHOD_CORR;
 819   factor.missing_type = MISS_LISTWISE;
 820   factor.exclude = MV_ANY;
 821   factor.print = PRINT_INITIAL | PRINT_EXTRACTION | PRINT_ROTATION;
 822   factor.extraction = EXTRACTION_PC;
 823   factor.n_factors = 0;
 824   factor.min_eigen = SYSMIS;
 825   factor.iterations = 25;
 826   factor.econverge = 0.001;
 827
 828   factor.blank = 0;
 829   factor.sort = false;
 830   factor.plot = 0;
 831   factor.rotation = ROT_VARIMAX;
 832
 833   factor.rconverge = 0.0001;
 834
 835   factor.wv = dict_get_weight (dict);
 836
 837   lex_match (lexer, T_SLASH);
 838
 839   if (!lex_force_match_id (lexer, "VARIABLES"))
 840     {
 841       goto error;
 842     }
 843
 844   lex_match (lexer, T_EQUALS);
 845
 846   if (!parse_variables_const (lexer, dict, &factor.vars, &factor.n_vars,
 847                               PV_NO_DUPLICATE | PV_NUMERIC))
 848     goto error;
 849
 850   if (factor.n_vars < 2)
 851     msg (MW, _("Factor analysis on a single variable is not useful."));
 852
 853   while (lex_token (lexer) != T_ENDCMD)
 854     {
 855       lex_match (lexer, T_SLASH);
 856
 857       if (lex_match_id (lexer, "PLOT"))
 858         {
 859           lex_match (lexer, T_EQUALS);
 860           while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
 861             {
 862               if (lex_match_id (lexer, "EIGEN"))
 863                 {
 864                   factor.plot |= PLOT_SCREE;
 865                 }
 866 #if FACTOR_FULLY_IMPLEMENTED
 867               else if (lex_match_id (lexer, "ROTATION"))
 868                 {
 869                 }
 870 #endif
 871               else
 872                 {
 873                   lex_error (lexer, NULL);
 874                   goto error;
 875                 }
 876             }
 877         }
 878       else if (lex_match_id (lexer, "METHOD"))
 879         {
 880           lex_match (lexer, T_EQUALS);
 881           while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
 882             {
 883               if (lex_match_id (lexer, "COVARIANCE"))
 884                 {
 885                   factor.method = METHOD_COV;
 886                 }
 887               else if (lex_match_id (lexer, "CORRELATION"))
 888                 {
 889                   factor.method = METHOD_CORR;
 890                 }
 891               else
 892                 {
 893                   lex_error (lexer, NULL);
 894                   goto error;
 895                 }
 896             }
 897         }
 898       else if (lex_match_id (lexer, "ROTATION"))
 899         {
 900           lex_match (lexer, T_EQUALS);
 901           while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
 902             {
 903               /* VARIMAX and DEFAULT are defaults */
 904               if (lex_match_id (lexer, "VARIMAX") || lex_match_id (lexer, "DEFAULT"))
 905                 {
 906                   factor.rotation = ROT_VARIMAX;
 907                 }
 908               else if (lex_match_id (lexer, "EQUAMAX"))
 909                 {
 910                   factor.rotation = ROT_EQUAMAX;
 911                 }
 912               else if (lex_match_id (lexer, "QUARTIMAX"))
 913                 {
 914                   factor.rotation = ROT_QUARTIMAX;
 915                 }
 916               else if (lex_match_id (lexer, "NOROTATE"))
 917                 {
 918                   factor.rotation = ROT_NONE;
 919                 }
 920               else
 921                 {
 922                   lex_error (lexer, NULL);
 923                   goto error;
 924                 }
 925             }
 926         }
 927       else if (lex_match_id (lexer, "CRITERIA"))
 928         {
 929           lex_match (lexer, T_EQUALS);
 930           while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
 931             {
 932               if (lex_match_id (lexer, "FACTORS"))
 933                 {
 934                   if ( lex_force_match (lexer, T_LPAREN))
 935                     {
 936                       lex_force_int (lexer);
 937                       factor.n_factors = lex_integer (lexer);
 938                       lex_get (lexer);
 939                       lex_force_match (lexer, T_RPAREN);
 940                     }
 941                 }
 942               else if (lex_match_id (lexer, "MINEIGEN"))
 943                 {
 944                   if ( lex_force_match (lexer, T_LPAREN))
 945                     {
 946                       lex_force_num (lexer);
 947                       factor.min_eigen = lex_number (lexer);
 948                       lex_get (lexer);
 949                       lex_force_match (lexer, T_RPAREN);
 950                     }
 951                 }
 952               else if (lex_match_id (lexer, "ECONVERGE"))
 953                 {
 954                   if ( lex_force_match (lexer, T_LPAREN))
 955                     {
 956                       lex_force_num (lexer);
 957                       factor.econverge = lex_number (lexer);
 958                       lex_get (lexer);
 959                       lex_force_match (lexer, T_RPAREN);
 960                     }
 961                 }
 962               else if (lex_match_id (lexer, "RCONVERGE"))
 963                 {
 964                   if ( lex_force_match (lexer, T_LPAREN))
 965                     {
 966                       lex_force_num (lexer);
 967                       factor.rconverge = lex_number (lexer);
 968                       lex_get (lexer);
 969                       lex_force_match (lexer, T_RPAREN);
 970                     }
 971                 }
 972               else if (lex_match_id (lexer, "ITERATE"))
 973                 {
 974                   if ( lex_force_match (lexer, T_LPAREN))
 975                     {
 976                       lex_force_int (lexer);
 977                       factor.iterations = lex_integer (lexer);
 978                       lex_get (lexer);
 979                       lex_force_match (lexer, T_RPAREN);
 980                     }
 981                 }
 982               else if (lex_match_id (lexer, "DEFAULT"))
 983                 {
 984                   factor.n_factors = 0;
 985                   factor.min_eigen = 1;
 986                   factor.iterations = 25;
 987                 }
 988               else
 989                 {
 990                   lex_error (lexer, NULL);
 991                   goto error;
 992                 }
 993             }
 994         }
 995       else if (lex_match_id (lexer, "EXTRACTION"))
 996         {
 997           extraction_seen = true;
 998           lex_match (lexer, T_EQUALS);
 999           while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
1000             {
1001               if (lex_match_id (lexer, "PAF"))
1002                 {
1003                   factor.extraction = EXTRACTION_PAF;
1004                 }
1005               else if (lex_match_id (lexer, "PC"))
1006                 {
1007                   factor.extraction = EXTRACTION_PC;
1008                 }
1009               else if (lex_match_id (lexer, "PA1"))
1010                 {
1011                   factor.extraction = EXTRACTION_PC;
1012                 }
1013               else if (lex_match_id (lexer, "DEFAULT"))
1014                 {
1015                   factor.extraction = EXTRACTION_PC;
1016                 }
1017               else
1018                 {
1019                   lex_error (lexer, NULL);
1020                   goto error;
1021                 }
1022             }
1023         }
1024       else if (lex_match_id (lexer, "FORMAT"))
1025         {
1026           lex_match (lexer, T_EQUALS);
1027           while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
1028             {
1029               if (lex_match_id (lexer, "SORT"))
1030                 {
1031                   factor.sort = true;
1032                 }
1033               else if (lex_match_id (lexer, "BLANK"))
1034                 {
1035                   if ( lex_force_match (lexer, T_LPAREN))
1036                     {
1037                       lex_force_num (lexer);
1038                       factor.blank = lex_number (lexer);
1039                       lex_get (lexer);
1040                       lex_force_match (lexer, T_RPAREN);
1041                     }
1042                 }
1043               else if (lex_match_id (lexer, "DEFAULT"))
1044                 {
1045                   factor.blank = 0;
1046                   factor.sort = false;
1047                 }
1048               else
1049                 {
1050                   lex_error (lexer, NULL);
1051                   goto error;
1052                 }
1053             }
1054         }
1055       else if (lex_match_id (lexer, "PRINT"))
1056         {
1057           factor.print = 0;
1058           lex_match (lexer, T_EQUALS);
1059           while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
1060             {
1061               if (lex_match_id (lexer, "UNIVARIATE"))
1062                 {
1063                   factor.print |= PRINT_UNIVARIATE;
1064                 }
1065               else if (lex_match_id (lexer, "DET"))
1066                 {
1067                   factor.print |= PRINT_DETERMINANT;
1068                 }
1069 #if FACTOR_FULLY_IMPLEMENTED
1070               else if (lex_match_id (lexer, "INV"))
1071                 {
1072                 }
1073               else if (lex_match_id (lexer, "AIC"))
1074                 {
1075                 }
1076 #endif
1077               else if (lex_match_id (lexer, "SIG"))
1078                 {
1079                   factor.print |= PRINT_SIG;
1080                 }
1081               else if (lex_match_id (lexer, "CORRELATION"))
1082                 {
1083                   factor.print |= PRINT_CORRELATION;
1084                 }
1085 #if FACTOR_FULLY_IMPLEMENTED
1086               else if (lex_match_id (lexer, "COVARIANCE"))
1087                 {
1088                 }
1089 #endif
1090               else if (lex_match_id (lexer, "ROTATION"))
1091                 {
1092                   factor.print |= PRINT_ROTATION;
1093                 }
1094               else if (lex_match_id (lexer, "EXTRACTION"))
1095                 {
1096                   factor.print |= PRINT_EXTRACTION;
1097                 }
1098               else if (lex_match_id (lexer, "INITIAL"))
1099                 {
1100                   factor.print |= PRINT_INITIAL;
1101                 }
1102               else if (lex_match_id (lexer, "KMO"))
1103                 {
1104                   factor.print |= PRINT_KMO;
1105                 }
1106 #if FACTOR_FULLY_IMPLEMENTED
1107               else if (lex_match_id (lexer, "REPR"))
1108                 {
1109                 }
1110               else if (lex_match_id (lexer, "FSCORE"))
1111                 {
1112                 }
1113 #endif
1114               else if (lex_match (lexer, T_ALL))
1115                 {
1116                   factor.print = 0xFFFF;
1117                 }
1118               else if (lex_match_id (lexer, "DEFAULT"))
1119                 {
1120                   factor.print |= PRINT_INITIAL ;
1121                   factor.print |= PRINT_EXTRACTION ;
1122                   factor.print |= PRINT_ROTATION ;
1123                 }
1124               else
1125                 {
1126                   lex_error (lexer, NULL);
1127                   goto error;
1128                 }
1129             }
1130         }
1131       else if (lex_match_id (lexer, "MISSING"))
1132         {
1133           lex_match (lexer, T_EQUALS);
1134           while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
1135             {
1136               if (lex_match_id (lexer, "INCLUDE"))
1137                 {
1138                   factor.exclude = MV_SYSTEM;
1139                 }
1140               else if (lex_match_id (lexer, "EXCLUDE"))
1141                 {
1142                   factor.exclude = MV_ANY;
1143                 }
1144               else if (lex_match_id (lexer, "LISTWISE"))
1145                 {
1146                   factor.missing_type = MISS_LISTWISE;
1147                 }
1148               else if (lex_match_id (lexer, "PAIRWISE"))
1149                 {
1150                   factor.missing_type = MISS_PAIRWISE;
1151                 }
1152               else if (lex_match_id (lexer, "MEANSUB"))
1153                 {
1154                   factor.missing_type = MISS_MEANSUB;
1155                 }
1156               else
1157                 {
1158                   lex_error (lexer, NULL);
1159                   goto error;
1160                 }
1161             }
1162         }
1163       else
1164         {
1165           lex_error (lexer, NULL);
1166           goto error;
1167         }
1168     }
1169
1170   if ( factor.rotation == ROT_NONE )
1171     factor.print &= ~PRINT_ROTATION;
1172
1173   if ( ! run_factor (ds, &factor))
1174     goto error;
1175
1176   free (factor.vars);
1177   return CMD_SUCCESS;
1178
1179  error:
1180   free (factor.vars);
1181   return CMD_FAILURE;
1182 }
1183
1184 static void do_factor (const struct cmd_factor *factor, struct casereader *group);
1185
1186
1187 static bool
1188 run_factor (struct dataset *ds, const struct cmd_factor *factor)
1189 {
1190   struct dictionary *dict = dataset_dict (ds);
1191   bool ok;
1192   struct casereader *group;
1193
1194   struct casegrouper *grouper = casegrouper_create_splits (proc_open (ds), dict);
1195
1196   while (casegrouper_get_next_group (grouper, &group))
1197     {
1198       if ( factor->missing_type == MISS_LISTWISE )
1199         group  = casereader_create_filter_missing (group, factor->vars, factor->n_vars,
1200                                                    factor->exclude,
1201                                                    NULL,  NULL);
1202       do_factor (factor, group);
1203     }
1204
1205   ok = casegrouper_destroy (grouper);
1206   ok = proc_commit (ds) && ok;
1207
1208   return ok;
1209 }
1210
1211
1212 /* Return the communality of variable N, calculated to N_FACTORS */
1213 static double
1214 the_communality (const gsl_matrix *evec, const gsl_vector *eval, int n, int n_factors)
1215 {
1216   size_t i;
1217
1218   double comm = 0;
1219
1220   assert (n >= 0);
1221   assert (n < eval->size);
1222   assert (n < evec->size1);
1223   assert (n_factors <= eval->size);
1224
1225   for (i = 0 ; i < n_factors; ++i)
1226     {
1227       double evali = fabs (gsl_vector_get (eval, i));
1228
1229       double eveci = gsl_matrix_get (evec, n, i);
1230
1231       comm += pow2 (eveci) * evali;
1232     }
1233
1234   return comm;
1235 }
1236
1237 /* Return the communality of variable N, calculated to N_FACTORS */
1238 static double
1239 communality (struct idata *idata, int n, int n_factors)
1240 {
1241   return the_communality (idata->evec, idata->eval, n, n_factors);
1242 }
1243
1244
1245 static void
1246 show_scree (const struct cmd_factor *f, struct idata *idata)
1247 {
1248   struct scree *s;
1249   const char *label ;
1250
1251   if ( !(f->plot & PLOT_SCREE) )
1252     return;
1253
1254
1255   label = f->extraction == EXTRACTION_PC ? _("Component Number") : _("Factor Number");
1256
1257   s = scree_create (idata->eval, label);
1258
1259   scree_submit (s);
1260 }
1261
1262 static void
1263 show_communalities (const struct cmd_factor * factor,
1264                     const gsl_vector *initial, const gsl_vector *extracted)
1265 {
1266   int i;
1267   int c = 0;
1268   const int heading_columns = 1;
1269   int nc = heading_columns;
1270   const int heading_rows = 1;
1271   const int nr = heading_rows + factor->n_vars;
1272   struct tab_table *t;
1273
1274   if (factor->print & PRINT_EXTRACTION)
1275     nc++;
1276
1277   if (factor->print & PRINT_INITIAL)
1278     nc++;
1279
1280   /* No point having a table with only headings */
1281   if (nc <= 1)
1282     return;
1283
1284   t = tab_create (nc, nr);
1285
1286   tab_title (t, _("Communalities"));
1287
1288   tab_headers (t, heading_columns, 0, heading_rows, 0);
1289
1290   c = 1;
1291   if (factor->print & PRINT_INITIAL)
1292     tab_text (t, c++, 0, TAB_CENTER | TAT_TITLE, _("Initial"));
1293
1294   if (factor->print & PRINT_EXTRACTION)
1295     tab_text (t, c++, 0, TAB_CENTER | TAT_TITLE, _("Extraction"));
1296
1297   /* Outline the box */
1298   tab_box (t,
1299            TAL_2, TAL_2,
1300            -1, -1,
1301            0, 0,
1302            nc - 1, nr - 1);
1303
1304   /* Vertical lines */
1305   tab_box (t,
1306            -1, -1,
1307            -1, TAL_1,
1308            heading_columns, 0,
1309            nc - 1, nr - 1);
1310
1311   tab_hline (t, TAL_1, 0, nc - 1, heading_rows);
1312   tab_vline (t, TAL_2, heading_columns, 0, nr - 1);
1313
1314   for (i = 0 ; i < factor->n_vars; ++i)
1315     {
1316       c = 0;
1317       tab_text (t, c++, i + heading_rows, TAT_TITLE, var_to_string (factor->vars[i]));
1318
1319       if (factor->print & PRINT_INITIAL)
1320         tab_double (t, c++, i + heading_rows, 0, gsl_vector_get (initial, i), NULL);
1321
1322       if (factor->print & PRINT_EXTRACTION)
1323         tab_double (t, c++, i + heading_rows, 0, gsl_vector_get (extracted, i), NULL);
1324     }
1325
1326   tab_submit (t);
1327 }
1328
1329
1330 static void
1331 show_factor_matrix (const struct cmd_factor *factor, struct idata *idata, const char *title, const gsl_matrix *fm)
1332 {
1333   int i;
1334   const int n_factors = idata->n_extractions;
1335
1336   const int heading_columns = 1;
1337   const int heading_rows = 2;
1338   const int nr = heading_rows + factor->n_vars;
1339   const int nc = heading_columns + n_factors;
1340   gsl_permutation *perm;
1341
1342   struct tab_table *t = tab_create (nc, nr);
1343
1344   /*
1345   if ( factor->extraction == EXTRACTION_PC )
1346     tab_title (t, _("Component Matrix"));
1347   else
1348     tab_title (t, _("Factor Matrix"));
1349   */
1350
1351   tab_title (t, "%s", title);
1352
1353   tab_headers (t, heading_columns, 0, heading_rows, 0);
1354
1355   if ( factor->extraction == EXTRACTION_PC )
1356     tab_joint_text (t,
1357                     1, 0,
1358                     nc - 1, 0,
1359                     TAB_CENTER | TAT_TITLE, _("Component"));
1360   else
1361     tab_joint_text (t,
1362                     1, 0,
1363                     nc - 1, 0,
1364                     TAB_CENTER | TAT_TITLE, _("Factor"));
1365
1366
1367   tab_hline (t, TAL_1, heading_columns, nc - 1, 1);
1368
1369
1370   /* Outline the box */
1371   tab_box (t,
1372            TAL_2, TAL_2,
1373            -1, -1,
1374            0, 0,
1375            nc - 1, nr - 1);
1376
1377   /* Vertical lines */
1378   tab_box (t,
1379            -1, -1,
1380            -1, TAL_1,
1381            heading_columns, 1,
1382            nc - 1, nr - 1);
1383
1384   tab_hline (t, TAL_1, 0, nc - 1, heading_rows);
1385   tab_vline (t, TAL_2, heading_columns, 0, nr - 1);
1386
1387
1388   /* Initialise to the identity permutation */
1389   perm = gsl_permutation_calloc (factor->n_vars);
1390
1391   if ( factor->sort)
1392     sort_matrix_indirect (fm, perm);
1393
1394   for (i = 0 ; i < n_factors; ++i)
1395     {
1396       tab_text_format (t, heading_columns + i, 1, TAB_CENTER | TAT_TITLE, _("%d"), i + 1);
1397     }
1398
1399   for (i = 0 ; i < factor->n_vars; ++i)
1400     {
1401       int j;
1402       const int matrix_row = perm->data[i];
1403       tab_text (t, 0, i + heading_rows, TAT_TITLE, var_to_string (factor->vars[matrix_row]));
1404
1405       for (j = 0 ; j < n_factors; ++j)
1406         {
1407           double x = gsl_matrix_get (fm, matrix_row, j);
1408
1409           if ( fabs (x) < factor->blank)
1410             continue;
1411
1412           tab_double (t, heading_columns + j, heading_rows + i, 0, x, NULL);
1413         }
1414     }
1415
1416   gsl_permutation_free (perm);
1417
1418   tab_submit (t);
1419 }
1420
1421
1422 static void
1423 show_explained_variance (const struct cmd_factor * factor, struct idata *idata,
1424                          const gsl_vector *initial_eigenvalues,
1425                          const gsl_vector *extracted_eigenvalues,
1426                          const gsl_vector *rotated_loadings)
1427 {
1428   size_t i;
1429   int c = 0;
1430   const int heading_columns = 1;
1431   const int heading_rows = 2;
1432   const int nr = heading_rows + factor->n_vars;
1433
1434   struct tab_table *t ;
1435
1436   double i_total = 0.0;
1437   double i_cum = 0.0;
1438
1439   double e_total = 0.0;
1440   double e_cum = 0.0;
1441
1442   double r_cum = 0.0;
1443
1444   int nc = heading_columns;
1445
1446   if (factor->print & PRINT_EXTRACTION)
1447     nc += 3;
1448
1449   if (factor->print & PRINT_INITIAL)
1450     nc += 3;
1451
1452   if (factor->print & PRINT_ROTATION)
1453     nc += 3;
1454
1455   /* No point having a table with only headings */
1456   if ( nc <= heading_columns)
1457     return;
1458
1459   t = tab_create (nc, nr);
1460
1461   tab_title (t, _("Total Variance Explained"));
1462
1463   tab_headers (t, heading_columns, 0, heading_rows, 0);
1464
1465   /* Outline the box */
1466   tab_box (t,
1467            TAL_2, TAL_2,
1468            -1, -1,
1469            0, 0,
1470            nc - 1, nr - 1);
1471
1472   /* Vertical lines */
1473   tab_box (t,
1474            -1, -1,
1475            -1, TAL_1,
1476            heading_columns, 0,
1477            nc - 1, nr - 1);
1478
1479   tab_hline (t, TAL_1, 0, nc - 1, heading_rows);
1480   tab_hline (t, TAL_1, 1, nc - 1, 1);
1481
1482   tab_vline (t, TAL_2, heading_columns, 0, nr - 1);
1483
1484
1485   if ( factor->extraction == EXTRACTION_PC)
1486     tab_text (t, 0, 1, TAB_LEFT | TAT_TITLE, _("Component"));
1487   else
1488     tab_text (t, 0, 1, TAB_LEFT | TAT_TITLE, _("Factor"));
1489
1490   c = 1;
1491   if (factor->print & PRINT_INITIAL)
1492     {
1493       tab_joint_text (t, c, 0, c + 2, 0, TAB_CENTER | TAT_TITLE, _("Initial Eigenvalues"));
1494       c += 3;
1495     }
1496
1497   if (factor->print & PRINT_EXTRACTION)
1498     {
1499       tab_joint_text (t, c, 0, c + 2, 0, TAB_CENTER | TAT_TITLE, _("Extraction Sums of Squared Loadings"));
1500       c += 3;
1501     }
1502
1503   if (factor->print & PRINT_ROTATION)
1504     {
1505       tab_joint_text (t, c, 0, c + 2, 0, TAB_CENTER | TAT_TITLE, _("Rotation Sums of Squared Loadings"));
1506       c += 3;
1507     }
1508
1509   for (i = 0; i < (nc - heading_columns) / 3 ; ++i)
1510     {
1511       tab_text (t, i * 3 + 1, 1, TAB_CENTER | TAT_TITLE, _("Total"));
1512       /* xgettext:no-c-format */
1513       tab_text (t, i * 3 + 2, 1, TAB_CENTER | TAT_TITLE, _("% of Variance"));
1514       tab_text (t, i * 3 + 3, 1, TAB_CENTER | TAT_TITLE, _("Cumulative %"));
1515
1516       tab_vline (t, TAL_2, heading_columns + i * 3, 0, nr - 1);
1517     }
1518
1519   for (i = 0 ; i < initial_eigenvalues->size; ++i)
1520     i_total += gsl_vector_get (initial_eigenvalues, i);
1521
1522   if ( factor->extraction == EXTRACTION_PAF)
1523     {
1524       e_total = factor->n_vars;
1525     }
1526   else
1527     {
1528       e_total = i_total;
1529     }
1530
1531   for (i = 0 ; i < factor->n_vars; ++i)
1532     {
1533       const double i_lambda = gsl_vector_get (initial_eigenvalues, i);
1534       double i_percent = 100.0 * i_lambda / i_total ;
1535
1536       const double e_lambda = gsl_vector_get (extracted_eigenvalues, i);
1537       double e_percent = 100.0 * e_lambda / e_total ;
1538
1539       const double r_lambda = gsl_vector_get (rotated_loadings, i);
1540       double r_percent = 100.0 * r_lambda / e_total ;
1541
1542       c = 0;
1543
1544       tab_text_format (t, c++, i + heading_rows, TAB_LEFT | TAT_TITLE, _("%zu"), i + 1);
1545
1546       i_cum += i_percent;
1547       e_cum += e_percent;
1548       r_cum += r_percent;
1549
1550       /* Initial Eigenvalues */
1551       if (factor->print & PRINT_INITIAL)
1552       {
1553         tab_double (t, c++, i + heading_rows, 0, i_lambda, NULL);
1554         tab_double (t, c++, i + heading_rows, 0, i_percent, NULL);
1555         tab_double (t, c++, i + heading_rows, 0, i_cum, NULL);
1556       }
1557
1558
1559       if (factor->print & PRINT_EXTRACTION)
1560         {
1561           if (i < idata->n_extractions)
1562             {
1563               /* Sums of squared loadings */
1564               tab_double (t, c++, i + heading_rows, 0, e_lambda, NULL);
1565               tab_double (t, c++, i + heading_rows, 0, e_percent, NULL);
1566               tab_double (t, c++, i + heading_rows, 0, e_cum, NULL);
1567             }
1568         }
1569
1570       if (factor->print & PRINT_ROTATION)
1571         {
1572           if (i < idata->n_extractions)
1573             {
1574               tab_double (t, c++, i + heading_rows, 0, r_lambda, NULL);
1575               tab_double (t, c++, i + heading_rows, 0, r_percent, NULL);
1576               tab_double (t, c++, i + heading_rows, 0, r_cum, NULL);
1577             }
1578       }
1579
1580     }
1581
1582   tab_submit (t);
1583 }
1584
1585
1586 static void
1587 show_correlation_matrix (const struct cmd_factor *factor, const struct idata *idata)
1588 {
1589   struct tab_table *t ;
1590   size_t i, j;
1591   int y_pos_corr = -1;
1592   int y_pos_sig = -1;
1593   int suffix_rows = 0;
1594
1595   const int heading_rows = 1;
1596   const int heading_columns = 2;
1597
1598   int nc = heading_columns ;
1599   int nr = heading_rows ;
1600   int n_data_sets = 0;
1601
1602   if (factor->print & PRINT_CORRELATION)
1603     {
1604       y_pos_corr = n_data_sets;
1605       n_data_sets++;
1606       nc = heading_columns + factor->n_vars;
1607     }
1608
1609   if (factor->print & PRINT_SIG)
1610     {
1611       y_pos_sig = n_data_sets;
1612       n_data_sets++;
1613       nc = heading_columns + factor->n_vars;
1614     }
1615
1616   nr += n_data_sets * factor->n_vars;
1617
1618   if (factor->print & PRINT_DETERMINANT)
1619     suffix_rows = 1;
1620
1621   /* If the table would contain only headings, don't bother rendering it */
1622   if (nr <= heading_rows && suffix_rows == 0)
1623     return;
1624
1625   t = tab_create (nc, nr + suffix_rows);
1626
1627   tab_title (t, _("Correlation Matrix"));
1628
1629   tab_hline (t, TAL_1, 0, nc - 1, heading_rows);
1630
1631   if (nr > heading_rows)
1632     {
1633       tab_headers (t, heading_columns, 0, heading_rows, 0);
1634
1635       tab_vline (t, TAL_2, 2, 0, nr - 1);
1636
1637       /* Outline the box */
1638       tab_box (t,
1639                TAL_2, TAL_2,
1640                -1, -1,
1641                0, 0,
1642                nc - 1, nr - 1);
1643
1644       /* Vertical lines */
1645       tab_box (t,
1646                -1, -1,
1647                -1, TAL_1,
1648                heading_columns, 0,
1649                nc - 1, nr - 1);
1650
1651
1652       for (i = 0; i < factor->n_vars; ++i)
1653         tab_text (t, heading_columns + i, 0, TAT_TITLE, var_to_string (factor->vars[i]));
1654
1655
1656       for (i = 0 ; i < n_data_sets; ++i)
1657         {
1658           int y = heading_rows + i * factor->n_vars;
1659           size_t v;
1660           for (v = 0; v < factor->n_vars; ++v)
1661             tab_text (t, 1, y + v, TAT_TITLE, var_to_string (factor->vars[v]));
1662
1663           tab_hline (t, TAL_1, 0, nc - 1, y);
1664         }
1665
1666       if (factor->print & PRINT_CORRELATION)
1667         {
1668           const double y = heading_rows + y_pos_corr;
1669           tab_text (t, 0, y, TAT_TITLE, _("Correlations"));
1670
1671           for (i = 0; i < factor->n_vars; ++i)
1672             {
1673               for (j = 0; j < factor->n_vars; ++j)
1674                 tab_double (t, heading_columns + i,  y + j, 0, gsl_matrix_get (idata->corr, i, j), NULL);
1675             }
1676         }
1677
1678       if (factor->print & PRINT_SIG)
1679         {
1680           const double y = heading_rows + y_pos_sig * factor->n_vars;
1681           tab_text (t, 0, y, TAT_TITLE, _("Sig. (1-tailed)"));
1682
1683           for (i = 0; i < factor->n_vars; ++i)
1684             {
1685               for (j = 0; j < factor->n_vars; ++j)
1686                 {
1687                   double rho = gsl_matrix_get (idata->corr, i, j);
1688                   double w = gsl_matrix_get (idata->n, i, j);
1689
1690                   if (i == j)
1691                     continue;
1692
1693                   tab_double (t, heading_columns + i,  y + j, 0, significance_of_correlation (rho, w), NULL);
1694                 }
1695             }
1696         }
1697     }
1698
1699   if (factor->print & PRINT_DETERMINANT)
1700     {
1701       tab_text (t, 0, nr, TAB_LEFT | TAT_TITLE, _("Determinant"));
1702
1703       tab_double (t, 1, nr, 0, idata->detR, NULL);
1704     }
1705
1706   tab_submit (t);
1707 }
1708
1709
1710
1711 static void
1712 do_factor (const struct cmd_factor *factor, struct casereader *r)
1713 {
1714   struct ccase *c;
1715   const gsl_matrix *var_matrix;
1716   const gsl_matrix *mean_matrix;
1717
1718   const gsl_matrix *analysis_matrix;
1719   struct idata *idata = idata_alloc (factor->n_vars);
1720
1721   struct covariance *cov = covariance_1pass_create (factor->n_vars, factor->vars,
1722                                               factor->wv, factor->exclude);
1723
1724   for ( ; (c = casereader_read (r) ); case_unref (c))
1725     {
1726       covariance_accumulate (cov, c);
1727     }
1728
1729   idata->cov = covariance_calculate (cov);
1730
1731   if (idata->cov == NULL)
1732     {
1733       msg (MW, _("The dataset contains no complete observations. No analysis will be performed."));
1734       goto finish;
1735     }
1736
1737   var_matrix = covariance_moments (cov, MOMENT_VARIANCE);
1738   mean_matrix = covariance_moments (cov, MOMENT_MEAN);
1739   idata->n = covariance_moments (cov, MOMENT_NONE);
1740
1741   if ( factor->method == METHOD_CORR)
1742     {
1743       idata->corr = correlation_from_covariance (idata->cov, var_matrix);
1744
1745       analysis_matrix = idata->corr;
1746     }
1747   else
1748     analysis_matrix = idata->cov;
1749
1750   if (factor->print & PRINT_DETERMINANT
1751       || factor->print & PRINT_KMO)
1752     {
1753       int sign = 0;
1754
1755       const int size = idata->corr->size1;
1756       gsl_permutation *p = gsl_permutation_calloc (size);
1757       gsl_matrix *tmp = gsl_matrix_calloc (size, size);
1758       gsl_matrix_memcpy (tmp, idata->corr);
1759
1760       gsl_linalg_LU_decomp (tmp, p, &sign);
1761       idata->detR = gsl_linalg_LU_det (tmp, sign);
1762       gsl_permutation_free (p);
1763       gsl_matrix_free (tmp);
1764     }
1765
1766   if ( factor->print & PRINT_UNIVARIATE)
1767     {
1768       const struct fmt_spec *wfmt = factor->wv ? var_get_print_format (factor->wv) : & F_8_0;
1769       const int nc = 4;
1770       int i;
1771
1772       const int heading_columns = 1;
1773       const int heading_rows = 1;
1774
1775       const int nr = heading_rows + factor->n_vars;
1776
1777       struct tab_table *t = tab_create (nc, nr);
1778       tab_title (t, _("Descriptive Statistics"));
1779
1780       tab_headers (t, heading_columns, 0, heading_rows, 0);
1781
1782       /* Outline the box */
1783       tab_box (t,
1784                TAL_2, TAL_2,
1785                -1, -1,
1786                0, 0,
1787                nc - 1, nr - 1);
1788
1789       /* Vertical lines */
1790       tab_box (t,
1791                -1, -1,
1792                -1, TAL_1,
1793                heading_columns, 0,
1794                nc - 1, nr - 1);
1795
1796       tab_hline (t, TAL_1, 0, nc - 1, heading_rows);
1797       tab_vline (t, TAL_2, heading_columns, 0, nr - 1);
1798
1799       tab_text (t, 1, 0, TAB_CENTER | TAT_TITLE, _("Mean"));
1800       tab_text (t, 2, 0, TAB_CENTER | TAT_TITLE, _("Std. Deviation"));
1801       tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("Analysis N"));
1802
1803       for (i = 0 ; i < factor->n_vars; ++i)
1804         {
1805           const struct variable *v = factor->vars[i];
1806           tab_text (t, 0, i + heading_rows, TAB_LEFT | TAT_TITLE, var_to_string (v));
1807
1808           tab_double (t, 1, i + heading_rows, 0, gsl_matrix_get (mean_matrix, i, i), NULL);
1809           tab_double (t, 2, i + heading_rows, 0, sqrt (gsl_matrix_get (var_matrix, i, i)), NULL);
1810           tab_double (t, 3, i + heading_rows, 0, gsl_matrix_get (idata->n, i, i), wfmt);
1811         }
1812
1813       tab_submit (t);
1814     }
1815
1816   if (factor->print & PRINT_KMO)
1817     {
1818       int i;
1819       double sum_ssq_r = 0;
1820       double sum_ssq_a = 0;
1821
1822       double df = factor->n_vars * ( factor->n_vars - 1) / 2;
1823
1824       double w = 0;
1825
1826
1827       double xsq;
1828
1829       const int heading_columns = 2;
1830       const int heading_rows = 0;
1831
1832       const int nr = heading_rows + 4;
1833       const int nc = heading_columns + 1;
1834
1835       gsl_matrix *a, *x;
1836
1837       struct tab_table *t = tab_create (nc, nr);
1838       tab_title (t, _("KMO and Bartlett's Test"));
1839
1840       x  = clone_matrix (idata->corr);
1841       gsl_linalg_cholesky_decomp (x);
1842       gsl_linalg_cholesky_invert (x);
1843
1844       a = anti_image (x);
1845
1846       for (i = 0; i < x->size1; ++i)
1847         {
1848           sum_ssq_r += ssq_od_n (x, i);
1849           sum_ssq_a += ssq_od_n (a, i);
1850         }
1851
1852       gsl_matrix_free (a);
1853       gsl_matrix_free (x);
1854
1855       tab_headers (t, heading_columns, 0, heading_rows, 0);
1856
1857       /* Outline the box */
1858       tab_box (t,
1859                TAL_2, TAL_2,
1860                -1, -1,
1861                0, 0,
1862                nc - 1, nr - 1);
1863
1864       tab_vline (t, TAL_2, heading_columns, 0, nr - 1);
1865
1866       tab_text (t, 0, 0, TAT_TITLE | TAB_LEFT, _("Kaiser-Meyer-Olkin Measure of Sampling Adequacy"));
1867
1868       tab_double (t, 2, 0, 0, sum_ssq_r /  (sum_ssq_r + sum_ssq_a), NULL);
1869
1870       tab_text (t, 0, 1, TAT_TITLE | TAB_LEFT, _("Bartlett's Test of Sphericity"));
1871
1872       tab_text (t, 1, 1, TAT_TITLE, _("Approx. Chi-Square"));
1873       tab_text (t, 1, 2, TAT_TITLE, _("df"));
1874       tab_text (t, 1, 3, TAT_TITLE, _("Sig."));
1875
1876
1877       /* The literature doesn't say what to do for the value of W when
1878          missing values are involved.  The best thing I can think of
1879          is to take the mean average. */
1880       w = 0;
1881       for (i = 0; i < idata->n->size1; ++i)
1882         w += gsl_matrix_get (idata->n, i, i);
1883       w /= idata->n->size1;
1884
1885       xsq = w - 1 - (2 * factor->n_vars + 5) / 6.0;
1886       xsq *= -log (idata->detR);
1887
1888       tab_double (t, 2, 1, 0, xsq, NULL);
1889       tab_double (t, 2, 2, 0, df, &F_8_0);
1890       tab_double (t, 2, 3, 0, gsl_cdf_chisq_Q (xsq, df), NULL);
1891
1892
1893       tab_submit (t);
1894     }
1895
1896   show_correlation_matrix (factor, idata);
1897
1898   {
1899     gsl_eigen_symmv_workspace *workspace = gsl_eigen_symmv_alloc (factor->n_vars);
1900
1901     gsl_eigen_symmv (matrix_dup (analysis_matrix), idata->eval, idata->evec, workspace);
1902
1903     gsl_eigen_symmv_free (workspace);
1904   }
1905
1906   gsl_eigen_symmv_sort (idata->eval, idata->evec, GSL_EIGEN_SORT_ABS_DESC);
1907
1908   idata->n_extractions = n_extracted_factors (factor, idata);
1909
1910   if (idata->n_extractions == 0)
1911     {
1912       msg (MW, _("The FACTOR criteria result in zero factors extracted. Therefore no analysis will be performed."));
1913       goto finish;
1914     }
1915
1916   if (idata->n_extractions > factor->n_vars)
1917     {
1918       msg (MW, _("The FACTOR criteria result in more factors than variables, which is not meaningful. No analysis will be performed."));
1919       goto finish;
1920     }
1921
1922   {
1923     gsl_matrix *rotated_factors = NULL;
1924     gsl_vector *rotated_loadings = NULL;
1925
1926     const gsl_vector *extracted_eigenvalues = NULL;
1927     gsl_vector *initial_communalities = gsl_vector_alloc (factor->n_vars);
1928     gsl_vector *extracted_communalities = gsl_vector_alloc (factor->n_vars);
1929     size_t i;
1930     struct factor_matrix_workspace *fmw = factor_matrix_workspace_alloc (idata->msr->size, idata->n_extractions);
1931     gsl_matrix *factor_matrix = gsl_matrix_calloc (factor->n_vars, fmw->n_factors);
1932
1933     if ( factor->extraction == EXTRACTION_PAF)
1934       {
1935         gsl_vector *diff = gsl_vector_alloc (idata->msr->size);
1936         struct smr_workspace *ws = ws_create (analysis_matrix);
1937
1938         for (i = 0 ; i < factor->n_vars ; ++i)
1939           {
1940             double r2 = squared_multiple_correlation (analysis_matrix, i, ws);
1941
1942             gsl_vector_set (idata->msr, i, r2);
1943           }
1944         ws_destroy (ws);
1945
1946         gsl_vector_memcpy (initial_communalities, idata->msr);
1947
1948         for (i = 0; i < factor->iterations; ++i)
1949           {
1950             double min, max;
1951             gsl_vector_memcpy (diff, idata->msr);
1952
1953             iterate_factor_matrix (analysis_matrix, idata->msr, factor_matrix, fmw);
1954
1955             gsl_vector_sub (diff, idata->msr);
1956
1957             gsl_vector_minmax (diff, &min, &max);
1958
1959             if ( fabs (min) < factor->econverge && fabs (max) < factor->econverge)
1960               break;
1961           }
1962         gsl_vector_free (diff);
1963
1964
1965
1966         gsl_vector_memcpy (extracted_communalities, idata->msr);
1967         extracted_eigenvalues = fmw->eval;
1968       }
1969     else if (factor->extraction == EXTRACTION_PC)
1970       {
1971         for (i = 0; i < factor->n_vars; ++i)
1972           gsl_vector_set (initial_communalities, i, communality (idata, i, factor->n_vars));
1973
1974         gsl_vector_memcpy (extracted_communalities, initial_communalities);
1975
1976         iterate_factor_matrix (analysis_matrix, extracted_communalities, factor_matrix, fmw);
1977
1978
1979         extracted_eigenvalues = idata->eval;
1980       }
1981
1982
1983     show_communalities (factor, initial_communalities, extracted_communalities);
1984
1985
1986     if ( factor->rotation != ROT_NONE)
1987       {
1988         rotated_factors = gsl_matrix_calloc (factor_matrix->size1, factor_matrix->size2);
1989         rotated_loadings = gsl_vector_calloc (factor_matrix->size2);
1990
1991         rotate (factor, factor_matrix, extracted_communalities, rotated_factors, rotated_loadings);
1992       }
1993
1994     show_explained_variance (factor, idata, idata->eval, extracted_eigenvalues, rotated_loadings);
1995
1996     factor_matrix_workspace_free (fmw);
1997
1998     show_scree (factor, idata);
1999
2000     show_factor_matrix (factor, idata,
2001                         factor->extraction == EXTRACTION_PC ? _("Component Matrix") : _("Factor Matrix"),
2002                         factor_matrix);
2003
2004     if ( factor->rotation != ROT_NONE)
2005       {
2006         show_factor_matrix (factor, idata,
2007                             factor->extraction == EXTRACTION_PC ? _("Rotated Component Matrix") : _("Rotated Factor Matrix"),
2008                             rotated_factors);
2009
2010         gsl_matrix_free (rotated_factors);
2011       }
2012
2013
2014
2015     gsl_vector_free (initial_communalities);
2016     gsl_vector_free (extracted_communalities);
2017   }
2018
2019  finish:
2020
2021   idata_free (idata);
2022
2023   casereader_destroy (r);
2024 }
2025
2026
2027