X-Git-Url: https://pintos-os.org/cgi-bin/gitweb.cgi?a=blobdiff_plain;f=src%2Flanguage%2Fstats%2Ffactor.c;h=1145a48e66989064c18d48b4dfeb0812a538882a;hb=2be9bee9da6a2ce27715e58128569594319abfa2;hp=ca09aa5962e7f4fba8be9417ec6bd209f0412ad7;hpb=fbc7e9617f218331dbf6580c220f03a632a72a2e;p=pspp-builds.git diff --git a/src/language/stats/factor.c b/src/language/stats/factor.c index ca09aa59..1145a48e 100644 --- a/src/language/stats/factor.c +++ b/src/language/stats/factor.c @@ -1,5 +1,5 @@ /* PSPP - a program for statistical analysis. - Copyright (C) 2009 Free Software Foundation, Inc. + Copyright (C) 2009, 2010, 2011 Free Software Foundation, Inc. This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by @@ -16,7 +16,6 @@ #include - #include #include #include @@ -24,28 +23,25 @@ #include #include -#include - -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include -#include -#include - -#include -#include - -#include - +#include "data/casegrouper.h" +#include "data/casereader.h" +#include "data/casewriter.h" +#include "data/dataset.h" +#include "data/dictionary.h" +#include "data/format.h" +#include "data/subcase.h" +#include "language/command.h" +#include "language/lexer/lexer.h" +#include "language/lexer/value-parser.h" +#include "language/lexer/variable-parser.h" +#include "libpspp/message.h" +#include "libpspp/misc.h" +#include "math/correlation.h" +#include "math/covariance.h" +#include "math/moments.h" +#include "output/chart-item.h" +#include "output/charts/scree.h" +#include "output/tab.h" #include "gettext.h" #define _(msgid) gettext (msgid) @@ -70,6 +66,12 @@ enum extraction_method EXTRACTION_PAF, }; +enum plot_opts + { + PLOT_SCREE = 0x0001, + PLOT_ROTATION = 0x0002 + }; + enum print_opts { PRINT_UNIVARIATE = 0x0001, @@ -87,6 +89,52 @@ enum print_opts PRINT_FSCORE = 0x1000 }; +enum rotation_type + { + ROT_VARIMAX = 0, + ROT_EQUAMAX, + ROT_QUARTIMAX, + ROT_NONE + }; + +typedef void (*rotation_coefficients) (double *x, double *y, + double a, double b, double c, double d, + const gsl_matrix *loadings ); + + +static void +varimax_coefficients (double *x, double *y, + double a, double b, double c, double d, + const gsl_matrix *loadings ) +{ + *x = d - 2 * a * b / loadings->size1; + *y = c - (a * a - b * b) / loadings->size1; +} + +static void +equamax_coefficients (double *x, double *y, + double a, double b, double c, double d, + const gsl_matrix *loadings ) +{ + *x = d - loadings->size2 * a * b / loadings->size1; + *y = c - loadings->size2 * (a * a - b * b) / (2 * loadings->size1); +} + +static void +quartimax_coefficients (double *x, double *y, + double a UNUSED, double b UNUSED, double c, double d, + const gsl_matrix *loadings UNUSED) +{ + *x = d ; + *y = c ; +} + +static const rotation_coefficients rotation_coeff[3] = { + varimax_coefficients, + equamax_coefficients, + quartimax_coefficients +}; + struct cmd_factor { @@ -100,6 +148,8 @@ struct cmd_factor enum mv_class exclude; enum print_opts print; enum extraction_method extraction; + enum plot_opts plot; + enum rotation_type rotation; /* Extraction Criteria */ int n_factors; @@ -107,6 +157,8 @@ struct cmd_factor double econverge; int iterations; + double rconverge; + /* Format */ double blank; bool sort; @@ -116,6 +168,10 @@ struct idata { /* Intermediate values used in calculation */ + const gsl_matrix *corr ; /* The correlation matrix */ + gsl_matrix *cov ; /* The covariance matrix */ + const gsl_matrix *n ; /* Matrix of number of samples */ + gsl_vector *eval ; /* The eigenvalues */ gsl_matrix *evec ; /* The eigenvectors */ @@ -127,7 +183,7 @@ struct idata static struct idata * idata_alloc (size_t n_vars) { - struct idata *id = xmalloc (sizeof (*id)); + struct idata *id = xzalloc (sizeof (*id)); id->n_extractions = 0; id->msr = gsl_vector_alloc (n_vars); @@ -144,11 +200,14 @@ idata_free (struct idata *id) gsl_vector_free (id->msr); gsl_vector_free (id->eval); gsl_matrix_free (id->evec); + if (id->cov != NULL) + gsl_matrix_free (id->cov); free (id); } +#if 0 static void dump_matrix (const gsl_matrix *m) { @@ -187,6 +246,7 @@ dump_vector (const gsl_vector *v) } printf ("\n"); } +#endif static int @@ -280,7 +340,7 @@ ws_destroy (struct smr_workspace *ws) Return the square of the regression coefficient for VAR regressed against all other variables. */ static double -squared_multiple_correlation (const gsl_matrix *analysis_matrix, int var, struct smr_workspace *ws) +squared_multiple_correlation (const gsl_matrix *corr, int var, struct smr_workspace *ws) { /* For an explanation of what this is doing, see http://www.visualstatistics.net/Visual%20Statistics%20Multimedia/multiple_regression_analysis.htm @@ -289,7 +349,7 @@ squared_multiple_correlation (const gsl_matrix *analysis_matrix, int var, struct int signum = 0; gsl_matrix_view rxx; - gsl_matrix_memcpy (ws->m, analysis_matrix); + gsl_matrix_memcpy (ws->m, corr); gsl_matrix_swap_rows (ws->m, 0, var); gsl_matrix_swap_columns (ws->m, 0, var); @@ -449,13 +509,206 @@ sort_matrix_indirect (const gsl_matrix *input, gsl_permutation *perm) } +static void +drot_go (double phi, double *l0, double *l1) +{ + double r0 = cos (phi) * *l0 + sin (phi) * *l1; + double r1 = - sin (phi) * *l0 + cos (phi) * *l1; + + *l0 = r0; + *l1 = r1; +} + + +static gsl_matrix * +clone_matrix (const gsl_matrix *m) +{ + int j, k; + gsl_matrix *c = gsl_matrix_calloc (m->size1, m->size2); + + for (j = 0 ; j < c->size1; ++j) + { + for (k = 0 ; k < c->size2; ++k) + { + const double *v = gsl_matrix_const_ptr (m, j, k); + gsl_matrix_set (c, j, k, *v); + } + } + + return c; +} + + +static double +initial_sv (const gsl_matrix *fm) +{ + int j, k; + + double sv = 0.0; + for (j = 0 ; j < fm->size2; ++j) + { + double l4s = 0; + double l2s = 0; + + for (k = j + 1 ; k < fm->size2; ++k) + { + double lambda = gsl_matrix_get (fm, k, j); + double lambda_sq = lambda * lambda; + double lambda_4 = lambda_sq * lambda_sq; + + l4s += lambda_4; + l2s += lambda_sq; + } + sv += ( fm->size1 * l4s - (l2s * l2s) ) / (fm->size1 * fm->size1 ); + } + return sv; +} + +static void +rotate (const struct cmd_factor *cf, const gsl_matrix *unrot, + const gsl_vector *communalities, + gsl_matrix *result, + gsl_vector *rotated_loadings + ) +{ + int j, k; + int i; + double prev_sv; + + /* First get a normalised version of UNROT */ + gsl_matrix *normalised = gsl_matrix_calloc (unrot->size1, unrot->size2); + gsl_matrix *h_sqrt = gsl_matrix_calloc (communalities->size, communalities->size); + gsl_matrix *h_sqrt_inv ; + + /* H is the diagonal matrix containing the absolute values of the communalities */ + for (i = 0 ; i < communalities->size ; ++i) + { + double *ptr = gsl_matrix_ptr (h_sqrt, i, i); + *ptr = fabs (gsl_vector_get (communalities, i)); + } + + /* Take the square root of the communalities */ + gsl_linalg_cholesky_decomp (h_sqrt); + + + /* Save a copy of h_sqrt and invert it */ + h_sqrt_inv = clone_matrix (h_sqrt); + gsl_linalg_cholesky_decomp (h_sqrt_inv); + gsl_linalg_cholesky_invert (h_sqrt_inv); + + /* normalised vertion is H^{1/2} x UNROT */ + gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, h_sqrt_inv, unrot, 0.0, normalised); + + gsl_matrix_free (h_sqrt_inv); + + + /* Now perform the rotation iterations */ + + prev_sv = initial_sv (normalised); + for (i = 0 ; i < cf->iterations ; ++i) + { + double sv = 0.0; + for (j = 0 ; j < normalised->size2; ++j) + { + /* These variables relate to the convergence criterium */ + double l4s = 0; + double l2s = 0; + + for (k = j + 1 ; k < normalised->size2; ++k) + { + int p; + double a = 0.0; + double b = 0.0; + double c = 0.0; + double d = 0.0; + double x, y; + double phi; + + for (p = 0; p < normalised->size1; ++p) + { + double jv = gsl_matrix_get (normalised, p, j); + double kv = gsl_matrix_get (normalised, p, k); + + double u = jv * jv - kv * kv; + double v = 2 * jv * kv; + a += u; + b += v; + c += u * u - v * v; + d += 2 * u * v; + } + + rotation_coeff [cf->rotation] (&x, &y, a, b, c, d, normalised); + + phi = atan2 (x, y) / 4.0 ; + + /* Don't bother rotating if the angle is small */ + if ( fabs (sin (phi) ) <= pow (10.0, -15.0)) + continue; + + for (p = 0; p < normalised->size1; ++p) + { + double *lambda0 = gsl_matrix_ptr (normalised, p, j); + double *lambda1 = gsl_matrix_ptr (normalised, p, k); + drot_go (phi, lambda0, lambda1); + } + + /* Calculate the convergence criterium */ + { + double lambda = gsl_matrix_get (normalised, k, j); + double lambda_sq = lambda * lambda; + double lambda_4 = lambda_sq * lambda_sq; + + l4s += lambda_4; + l2s += lambda_sq; + } + } + sv += ( normalised->size1 * l4s - (l2s * l2s) ) / (normalised->size1 * normalised->size1 ); + } + + if ( fabs (sv - prev_sv) <= cf->rconverge) + break; + + prev_sv = sv; + } + + gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, + h_sqrt, normalised, 0.0, result); + + gsl_matrix_free (h_sqrt); + + + /* reflect negative sums and populate the rotated loadings vector*/ + for (i = 0 ; i < result->size2; ++i) + { + double ssq = 0.0; + double sum = 0.0; + for (j = 0 ; j < result->size1; ++j) + { + double s = gsl_matrix_get (result, j, i); + ssq += s * s; + sum += gsl_matrix_get (result, j, i); + } + + gsl_vector_set (rotated_loadings, i, ssq); + + if ( sum < 0 ) + for (j = 0 ; j < result->size1; ++j) + { + double *lambda = gsl_matrix_ptr (result, j, i); + *lambda = - *lambda; + } + } +} + + /* Get an approximation for the factor matrix into FACTORS, and the communalities into COMMUNALITIES. R is the matrix to be analysed. WS is a pointer to a structure which must have been initialised with factor_matrix_workspace_init. */ static void -iterate_factor_matrix (const gsl_matrix *r, gsl_vector *communalities, gsl_matrix *factors, struct factor_matrix_workspace *ws) +iterate_factor_matrix (const gsl_matrix *r, gsl_vector *communalities, gsl_matrix *factors, + struct factor_matrix_workspace *ws) { size_t i; gsl_matrix_view mv ; @@ -489,8 +742,7 @@ iterate_factor_matrix (const gsl_matrix *r, gsl_vector *communalities, gsl_matri /* Take the square root of gamma */ gsl_linalg_cholesky_decomp (ws->gamma); - gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, - 1.0, &mv.matrix, ws->gamma, 0.0, factors); + gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, &mv.matrix, ws->gamma, 0.0, factors); for (i = 0 ; i < r->size1 ; ++i) { @@ -511,6 +763,8 @@ cmd_factor (struct lexer *lexer, struct dataset *ds) const struct dictionary *dict = dataset_dict (ds); struct cmd_factor factor; + factor.n_vars = 0; + factor.vars = NULL; factor.method = METHOD_CORR; factor.missing_type = MISS_LISTWISE; factor.exclude = MV_ANY; @@ -520,40 +774,50 @@ cmd_factor (struct lexer *lexer, struct dataset *ds) factor.min_eigen = SYSMIS; factor.iterations = 25; factor.econverge = 0.001; + factor.blank = 0; factor.sort = false; + factor.plot = 0; + factor.rotation = ROT_VARIMAX; + + factor.rconverge = 0.0001; factor.wv = dict_get_weight (dict); - lex_match (lexer, '/'); + lex_match (lexer, T_SLASH); if (!lex_force_match_id (lexer, "VARIABLES")) { goto error; } - lex_match (lexer, '='); + lex_match (lexer, T_EQUALS); if (!parse_variables_const (lexer, dict, &factor.vars, &factor.n_vars, PV_NO_DUPLICATE | PV_NUMERIC)) goto error; - while (lex_token (lexer) != '.') + if (factor.n_vars < 2) + msg (MW, _("Factor analysis on a single variable is not useful.")); + + while (lex_token (lexer) != T_ENDCMD) { - lex_match (lexer, '/'); + lex_match (lexer, T_SLASH); -#if FACTOR_FULLY_IMPLEMENTED if (lex_match_id (lexer, "PLOT")) { - lex_match (lexer, '='); - while (lex_token (lexer) != '.' && lex_token (lexer) != '/') + lex_match (lexer, T_EQUALS); + while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH) { if (lex_match_id (lexer, "EIGEN")) { + factor.plot |= PLOT_SCREE; } +#if FACTOR_FULLY_IMPLEMENTED else if (lex_match_id (lexer, "ROTATION")) { } +#endif else { lex_error (lexer, NULL); @@ -561,12 +825,10 @@ cmd_factor (struct lexer *lexer, struct dataset *ds) } } } - else -#endif - if (lex_match_id (lexer, "METHOD")) + else if (lex_match_id (lexer, "METHOD")) { - lex_match (lexer, '='); - while (lex_token (lexer) != '.' && lex_token (lexer) != '/') + lex_match (lexer, T_EQUALS); + while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH) { if (lex_match_id (lexer, "COVARIANCE")) { @@ -583,17 +845,27 @@ cmd_factor (struct lexer *lexer, struct dataset *ds) } } } -#if FACTOR_FULLY_IMPLEMENTED else if (lex_match_id (lexer, "ROTATION")) { - lex_match (lexer, '='); - while (lex_token (lexer) != '.' && lex_token (lexer) != '/') + lex_match (lexer, T_EQUALS); + while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH) { - if (lex_match_id (lexer, "VARIMAX")) + /* VARIMAX and DEFAULT are defaults */ + if (lex_match_id (lexer, "VARIMAX") || lex_match_id (lexer, "DEFAULT")) { + factor.rotation = ROT_VARIMAX; } - else if (lex_match_id (lexer, "DEFAULT")) + else if (lex_match_id (lexer, "EQUAMAX")) + { + factor.rotation = ROT_EQUAMAX; + } + else if (lex_match_id (lexer, "QUARTIMAX")) + { + factor.rotation = ROT_QUARTIMAX; + } + else if (lex_match_id (lexer, "NOROTATE")) { + factor.rotation = ROT_NONE; } else { @@ -602,50 +874,59 @@ cmd_factor (struct lexer *lexer, struct dataset *ds) } } } -#endif else if (lex_match_id (lexer, "CRITERIA")) { - lex_match (lexer, '='); - while (lex_token (lexer) != '.' && lex_token (lexer) != '/') + lex_match (lexer, T_EQUALS); + while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH) { if (lex_match_id (lexer, "FACTORS")) { - if ( lex_force_match (lexer, '(')) + if ( lex_force_match (lexer, T_LPAREN)) { lex_force_int (lexer); factor.n_factors = lex_integer (lexer); lex_get (lexer); - lex_force_match (lexer, ')'); + lex_force_match (lexer, T_RPAREN); } } else if (lex_match_id (lexer, "MINEIGEN")) { - if ( lex_force_match (lexer, '(')) + if ( lex_force_match (lexer, T_LPAREN)) { lex_force_num (lexer); factor.min_eigen = lex_number (lexer); lex_get (lexer); - lex_force_match (lexer, ')'); + lex_force_match (lexer, T_RPAREN); } } else if (lex_match_id (lexer, "ECONVERGE")) { - if ( lex_force_match (lexer, '(')) + if ( lex_force_match (lexer, T_LPAREN)) { lex_force_num (lexer); factor.econverge = lex_number (lexer); lex_get (lexer); - lex_force_match (lexer, ')'); + lex_force_match (lexer, T_RPAREN); + } + } + else if (lex_match_id (lexer, "RCONVERGE")) + { + if ( lex_force_match (lexer, T_LPAREN)) + { + lex_force_num (lexer); + factor.rconverge = lex_number (lexer); + lex_get (lexer); + lex_force_match (lexer, T_RPAREN); } } else if (lex_match_id (lexer, "ITERATE")) { - if ( lex_force_match (lexer, '(')) + if ( lex_force_match (lexer, T_LPAREN)) { lex_force_int (lexer); factor.iterations = lex_integer (lexer); lex_get (lexer); - lex_force_match (lexer, ')'); + lex_force_match (lexer, T_RPAREN); } } else if (lex_match_id (lexer, "DEFAULT")) @@ -664,8 +945,8 @@ cmd_factor (struct lexer *lexer, struct dataset *ds) else if (lex_match_id (lexer, "EXTRACTION")) { extraction_seen = true; - lex_match (lexer, '='); - while (lex_token (lexer) != '.' && lex_token (lexer) != '/') + lex_match (lexer, T_EQUALS); + while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH) { if (lex_match_id (lexer, "PAF")) { @@ -692,8 +973,8 @@ cmd_factor (struct lexer *lexer, struct dataset *ds) } else if (lex_match_id (lexer, "FORMAT")) { - lex_match (lexer, '='); - while (lex_token (lexer) != '.' && lex_token (lexer) != '/') + lex_match (lexer, T_EQUALS); + while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH) { if (lex_match_id (lexer, "SORT")) { @@ -701,12 +982,12 @@ cmd_factor (struct lexer *lexer, struct dataset *ds) } else if (lex_match_id (lexer, "BLANK")) { - if ( lex_force_match (lexer, '(')) + if ( lex_force_match (lexer, T_LPAREN)) { lex_force_num (lexer); factor.blank = lex_number (lexer); lex_get (lexer); - lex_force_match (lexer, ')'); + lex_force_match (lexer, T_RPAREN); } } else if (lex_match_id (lexer, "DEFAULT")) @@ -724,8 +1005,8 @@ cmd_factor (struct lexer *lexer, struct dataset *ds) else if (lex_match_id (lexer, "PRINT")) { factor.print = 0; - lex_match (lexer, '='); - while (lex_token (lexer) != '.' && lex_token (lexer) != '/') + lex_match (lexer, T_EQUALS); + while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH) { if (lex_match_id (lexer, "UNIVARIATE")) { @@ -742,13 +1023,17 @@ cmd_factor (struct lexer *lexer, struct dataset *ds) else if (lex_match_id (lexer, "AIC")) { } +#endif else if (lex_match_id (lexer, "SIG")) { + factor.print |= PRINT_SIG; } - else if (lex_match_id (lexer, "COVARIANCE")) + else if (lex_match_id (lexer, "CORRELATION")) { + factor.print |= PRINT_CORRELATION; } - else if (lex_match_id (lexer, "CORRELATION")) +#if FACTOR_FULLY_IMPLEMENTED + else if (lex_match_id (lexer, "COVARIANCE")) { } #endif @@ -794,8 +1079,8 @@ cmd_factor (struct lexer *lexer, struct dataset *ds) } else if (lex_match_id (lexer, "MISSING")) { - lex_match (lexer, '='); - while (lex_token (lexer) != '.' && lex_token (lexer) != '/') + lex_match (lexer, T_EQUALS); + while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH) { if (lex_match_id (lexer, "INCLUDE")) { @@ -831,6 +1116,9 @@ cmd_factor (struct lexer *lexer, struct dataset *ds) } } + if ( factor.rotation == ROT_NONE ) + factor.print &= ~PRINT_ROTATION; + if ( ! run_factor (ds, &factor)) goto error; @@ -903,6 +1191,22 @@ communality (struct idata *idata, int n, int n_factors) } +static void +show_scree (const struct cmd_factor *f, struct idata *idata) +{ + struct scree *s; + const char *label ; + + if ( !(f->plot & PLOT_SCREE) ) + return; + + + label = f->extraction == EXTRACTION_PC ? _("Component Number") : _("Factor Number"); + + s = scree_create (idata->eval, label); + + scree_submit (s); +} static void show_communalities (const struct cmd_factor * factor, @@ -926,12 +1230,10 @@ show_communalities (const struct cmd_factor * factor, if (nc <= 1) return; - t = tab_create (nc, nr, 0); + t = tab_create (nc, nr); tab_title (t, _("Communalities")); - tab_dim (t, tab_natural_dimensions, NULL); - tab_headers (t, heading_columns, 0, heading_rows, 0); c = 1; @@ -975,10 +1277,10 @@ show_communalities (const struct cmd_factor * factor, static void -show_factor_matrix (const struct cmd_factor *factor, struct idata *idata, const gsl_matrix *fm) +show_factor_matrix (const struct cmd_factor *factor, struct idata *idata, const char *title, const gsl_matrix *fm) { int i; - const int n_factors = n_extracted_factors (factor, idata); + const int n_factors = idata->n_extractions; const int heading_columns = 1; const int heading_rows = 2; @@ -986,14 +1288,16 @@ show_factor_matrix (const struct cmd_factor *factor, struct idata *idata, const const int nc = heading_columns + n_factors; gsl_permutation *perm; - struct tab_table *t = tab_create (nc, nr, 0); + struct tab_table *t = tab_create (nc, nr); + /* if ( factor->extraction == EXTRACTION_PC ) tab_title (t, _("Component Matrix")); else tab_title (t, _("Factor Matrix")); + */ - tab_dim (t, tab_natural_dimensions, NULL); + tab_title (t, title); tab_headers (t, heading_columns, 0, heading_rows, 0); @@ -1030,11 +1334,6 @@ show_factor_matrix (const struct cmd_factor *factor, struct idata *idata, const tab_vline (t, TAL_2, heading_columns, 0, nr - 1); - { - gsl_vector_const_view r1 = gsl_matrix_const_row (fm, 0); - gsl_vector_const_view r2 = gsl_matrix_const_row (fm, 1); - } - /* Initialise to the identity permutation */ perm = gsl_permutation_calloc (factor->n_vars); @@ -1072,7 +1371,8 @@ show_factor_matrix (const struct cmd_factor *factor, struct idata *idata, const static void show_explained_variance (const struct cmd_factor * factor, struct idata *idata, const gsl_vector *initial_eigenvalues, - const gsl_vector *extracted_eigenvalues) + const gsl_vector *extracted_eigenvalues, + const gsl_vector *rotated_loadings) { size_t i; int c = 0; @@ -1088,6 +1388,8 @@ show_explained_variance (const struct cmd_factor * factor, struct idata *idata, double e_total = 0.0; double e_cum = 0.0; + double r_cum = 0.0; + int nc = heading_columns; if (factor->print & PRINT_EXTRACTION) @@ -1103,12 +1405,10 @@ show_explained_variance (const struct cmd_factor * factor, struct idata *idata, if ( nc <= heading_columns) return; - t = tab_create (nc, nr, 0); + t = tab_create (nc, nr); tab_title (t, _("Total Variance Explained")); - tab_dim (t, tab_natural_dimensions, NULL); - tab_headers (t, heading_columns, 0, heading_rows, 0); /* Outline the box */ @@ -1158,6 +1458,7 @@ show_explained_variance (const struct cmd_factor * factor, struct idata *idata, for (i = 0; i < (nc - heading_columns) / 3 ; ++i) { tab_text (t, i * 3 + 1, 1, TAB_CENTER | TAT_TITLE, _("Total")); + /* xgettext:no-c-format */ tab_text (t, i * 3 + 2, 1, TAB_CENTER | TAT_TITLE, _("% of Variance")); tab_text (t, i * 3 + 3, 1, TAB_CENTER | TAT_TITLE, _("Cumulative %")); @@ -1176,7 +1477,6 @@ show_explained_variance (const struct cmd_factor * factor, struct idata *idata, e_total = i_total; } - for (i = 0 ; i < factor->n_vars; ++i) { const double i_lambda = gsl_vector_get (initial_eigenvalues, i); @@ -1185,12 +1485,16 @@ show_explained_variance (const struct cmd_factor * factor, struct idata *idata, const double e_lambda = gsl_vector_get (extracted_eigenvalues, i); double e_percent = 100.0 * e_lambda / e_total ; + const double r_lambda = gsl_vector_get (rotated_loadings, i); + double r_percent = 100.0 * r_lambda / e_total ; + c = 0; tab_text_format (t, c++, i + heading_rows, TAB_LEFT | TAT_TITLE, _("%d"), i + 1); i_cum += i_percent; e_cum += e_percent; + r_cum += r_percent; /* Initial Eigenvalues */ if (factor->print & PRINT_INITIAL) @@ -1200,9 +1504,10 @@ show_explained_variance (const struct cmd_factor * factor, struct idata *idata, tab_double (t, c++, i + heading_rows, 0, i_cum, NULL); } + if (factor->print & PRINT_EXTRACTION) { - if ( i < n_extracted_factors (factor, idata)) + if (i < idata->n_extractions) { /* Sums of squared loadings */ tab_double (t, c++, i + heading_rows, 0, e_lambda, NULL); @@ -1210,6 +1515,154 @@ show_explained_variance (const struct cmd_factor * factor, struct idata *idata, tab_double (t, c++, i + heading_rows, 0, e_cum, NULL); } } + + if (factor->print & PRINT_ROTATION) + { + if (i < idata->n_extractions) + { + tab_double (t, c++, i + heading_rows, 0, r_lambda, NULL); + tab_double (t, c++, i + heading_rows, 0, r_percent, NULL); + tab_double (t, c++, i + heading_rows, 0, r_cum, NULL); + } + } + + } + + tab_submit (t); +} + + +static void +show_correlation_matrix (const struct cmd_factor *factor, const struct idata *idata) +{ + struct tab_table *t ; + size_t i, j; + int y_pos_corr = -1; + int y_pos_sig = -1; + int suffix_rows = 0; + + const int heading_rows = 1; + const int heading_columns = 2; + + int nc = heading_columns ; + int nr = heading_rows ; + int n_data_sets = 0; + + if (factor->print & PRINT_CORRELATION) + { + y_pos_corr = n_data_sets; + n_data_sets++; + nc = heading_columns + factor->n_vars; + } + + if (factor->print & PRINT_SIG) + { + y_pos_sig = n_data_sets; + n_data_sets++; + nc = heading_columns + factor->n_vars; + } + + nr += n_data_sets * factor->n_vars; + + if (factor->print & PRINT_DETERMINANT) + suffix_rows = 1; + + /* If the table would contain only headings, don't bother rendering it */ + if (nr <= heading_rows && suffix_rows == 0) + return; + + t = tab_create (nc, nr + suffix_rows); + + tab_title (t, _("Correlation Matrix")); + + tab_hline (t, TAL_1, 0, nc - 1, heading_rows); + + if (nr > heading_rows) + { + tab_headers (t, heading_columns, 0, heading_rows, 0); + + tab_vline (t, TAL_2, 2, 0, nr - 1); + + /* Outline the box */ + tab_box (t, + TAL_2, TAL_2, + -1, -1, + 0, 0, + nc - 1, nr - 1); + + /* Vertical lines */ + tab_box (t, + -1, -1, + -1, TAL_1, + heading_columns, 0, + nc - 1, nr - 1); + + + for (i = 0; i < factor->n_vars; ++i) + tab_text (t, heading_columns + i, 0, TAT_TITLE, var_to_string (factor->vars[i])); + + + for (i = 0 ; i < n_data_sets; ++i) + { + int y = heading_rows + i * factor->n_vars; + size_t v; + for (v = 0; v < factor->n_vars; ++v) + tab_text (t, 1, y + v, TAT_TITLE, var_to_string (factor->vars[v])); + + tab_hline (t, TAL_1, 0, nc - 1, y); + } + + if (factor->print & PRINT_CORRELATION) + { + const double y = heading_rows + y_pos_corr; + tab_text (t, 0, y, TAT_TITLE, _("Correlations")); + + for (i = 0; i < factor->n_vars; ++i) + { + for (j = 0; j < factor->n_vars; ++j) + tab_double (t, heading_columns + i, y + j, 0, gsl_matrix_get (idata->corr, i, j), NULL); + } + } + + if (factor->print & PRINT_SIG) + { + const double y = heading_rows + y_pos_sig * factor->n_vars; + tab_text (t, 0, y, TAT_TITLE, _("Sig. (1-tailed)")); + + for (i = 0; i < factor->n_vars; ++i) + { + for (j = 0; j < factor->n_vars; ++j) + { + double rho = gsl_matrix_get (idata->corr, i, j); + double w = gsl_matrix_get (idata->n, i, j); + + if (i == j) + continue; + + tab_double (t, heading_columns + i, y + j, 0, significance_of_correlation (rho, w), NULL); + } + } + } + } + + if (factor->print & PRINT_DETERMINANT) + { + int sign = 0; + double det = 0.0; + + const int size = idata->corr->size1; + gsl_permutation *p = gsl_permutation_calloc (size); + gsl_matrix *tmp = gsl_matrix_calloc (size, size); + gsl_matrix_memcpy (tmp, idata->corr); + + gsl_linalg_LU_decomp (tmp, p, &sign); + det = gsl_linalg_LU_det (tmp, sign); + gsl_permutation_free (p); + gsl_matrix_free (tmp); + + + tab_text (t, 0, nr, TAB_LEFT | TAT_TITLE, _("Determinant")); + tab_double (t, 1, nr, 0, det, NULL); } tab_submit (t); @@ -1221,15 +1674,13 @@ static void do_factor (const struct cmd_factor *factor, struct casereader *r) { struct ccase *c; - const gsl_matrix *cov_matrix; const gsl_matrix *var_matrix; const gsl_matrix *mean_matrix; - const gsl_matrix *n_matrix; const gsl_matrix *analysis_matrix; - struct idata *idata; + struct idata *idata = idata_alloc (factor->n_vars); - struct covariance *cov = covariance_create (factor->n_vars, factor->vars, + struct covariance *cov = covariance_1pass_create (factor->n_vars, factor->vars, factor->wv, factor->exclude); for ( ; (c = casereader_read (r) ); case_unref (c)) @@ -1237,18 +1688,25 @@ do_factor (const struct cmd_factor *factor, struct casereader *r) covariance_accumulate (cov, c); } - cov_matrix = covariance_calculate (cov); + idata->cov = covariance_calculate (cov); + + if (idata->cov == NULL) + { + msg (MW, _("The dataset contains no complete observations. No analysis will be performed.")); + goto finish; + } var_matrix = covariance_moments (cov, MOMENT_VARIANCE); mean_matrix = covariance_moments (cov, MOMENT_MEAN); - n_matrix = covariance_moments (cov, MOMENT_NONE); + idata->n = covariance_moments (cov, MOMENT_NONE); if ( factor->method == METHOD_CORR) { - analysis_matrix = correlation_from_covariance (cov_matrix, var_matrix); + idata->corr = correlation_from_covariance (idata->cov, var_matrix); + analysis_matrix = idata->corr; } else - analysis_matrix = cov_matrix; + analysis_matrix = idata->cov; if ( factor->print & PRINT_UNIVARIATE) { @@ -1262,9 +1720,8 @@ do_factor (const struct cmd_factor *factor, struct casereader *r) const int nr = heading_rows + factor->n_vars; - struct tab_table *t = tab_create (nc, nr, 0); + struct tab_table *t = tab_create (nc, nr); tab_title (t, _("Descriptive Statistics")); - tab_dim (t, tab_natural_dimensions, NULL); tab_headers (t, heading_columns, 0, heading_rows, 0); @@ -1296,54 +1753,13 @@ do_factor (const struct cmd_factor *factor, struct casereader *r) tab_double (t, 1, i + heading_rows, 0, gsl_matrix_get (mean_matrix, i, i), NULL); tab_double (t, 2, i + heading_rows, 0, sqrt (gsl_matrix_get (var_matrix, i, i)), NULL); - tab_double (t, 3, i + heading_rows, 0, gsl_matrix_get (n_matrix, i, i), wfmt); + tab_double (t, 3, i + heading_rows, 0, gsl_matrix_get (idata->n, i, i), wfmt); } tab_submit (t); } - if ( factor->print & PRINT_DETERMINANT) - { - const int nc = 2; - const int heading_columns = 0; - const int heading_rows = 0; - const int nr = 1; - struct tab_table *t ; - - int sign = 0; - double det = 0.0; - const int size = analysis_matrix->size1; - gsl_permutation *p = gsl_permutation_calloc (size); - gsl_matrix *tmp = gsl_matrix_calloc (size, size); - - gsl_matrix_memcpy (tmp, analysis_matrix); - gsl_linalg_LU_decomp (tmp, p, &sign); - det = gsl_linalg_LU_det (tmp, sign); - gsl_permutation_free (p); - gsl_matrix_free (tmp); - - t = tab_create (nc, nr, 0); - - if ( factor->method == METHOD_CORR) - tab_title (t, _("Correlation Matrix")); - else - tab_title (t, _("Covariance Matrix")); - - tab_dim (t, tab_natural_dimensions, NULL); - - tab_headers (t, heading_columns, 0, heading_rows, 0); - - tab_hline (t, TAL_1, 0, nc - 1, heading_rows); - - tab_text (t, 0, 0, TAB_LEFT | TAT_TITLE, _("Determinant")); - tab_double (t, 1, 0, 0, det, NULL); - - tab_submit (t); - } - - - idata = idata_alloc (factor->n_vars); - + show_correlation_matrix (factor, idata); #if 1 { @@ -1357,12 +1773,29 @@ do_factor (const struct cmd_factor *factor, struct casereader *r) gsl_eigen_symmv_sort (idata->eval, idata->evec, GSL_EIGEN_SORT_ABS_DESC); #endif + idata->n_extractions = n_extracted_factors (factor, idata); + + if (idata->n_extractions == 0) + { + msg (MW, _("The FACTOR criteria result in zero factors extracted. Therefore no analysis will be performed.")); + goto finish; + } + + if (idata->n_extractions > factor->n_vars) + { + msg (MW, _("The FACTOR criteria result in more factors than variables, which is not meaningful. No analysis will be performed.")); + goto finish; + } + { + gsl_matrix *rotated_factors = NULL; + gsl_vector *rotated_loadings = NULL; + const gsl_vector *extracted_eigenvalues = NULL; gsl_vector *initial_communalities = gsl_vector_alloc (factor->n_vars); gsl_vector *extracted_communalities = gsl_vector_alloc (factor->n_vars); size_t i; - struct factor_matrix_workspace *fmw = factor_matrix_workspace_alloc (idata->msr->size, n_extracted_factors (factor, idata)); + struct factor_matrix_workspace *fmw = factor_matrix_workspace_alloc (idata->msr->size, idata->n_extractions); gsl_matrix *factor_matrix = gsl_matrix_calloc (factor->n_vars, fmw->n_factors); if ( factor->extraction == EXTRACTION_PAF) @@ -1396,34 +1829,67 @@ do_factor (const struct cmd_factor *factor, struct casereader *r) } gsl_vector_free (diff); + + gsl_vector_memcpy (extracted_communalities, idata->msr); extracted_eigenvalues = fmw->eval; } else if (factor->extraction == EXTRACTION_PC) { - for (i = 0 ; i < factor->n_vars; ++i) - { - gsl_vector_set (initial_communalities, i, communality (idata, i, factor->n_vars)); - } + for (i = 0; i < factor->n_vars; ++i) + gsl_vector_set (initial_communalities, i, communality (idata, i, factor->n_vars)); + gsl_vector_memcpy (extracted_communalities, initial_communalities); iterate_factor_matrix (analysis_matrix, extracted_communalities, factor_matrix, fmw); + + extracted_eigenvalues = idata->eval; } + show_communalities (factor, initial_communalities, extracted_communalities); - show_explained_variance (factor, idata, idata->eval, extracted_eigenvalues); + + if ( factor->rotation != ROT_NONE) + { + rotated_factors = gsl_matrix_calloc (factor_matrix->size1, factor_matrix->size2); + rotated_loadings = gsl_vector_calloc (factor_matrix->size2); + + rotate (factor, factor_matrix, extracted_communalities, rotated_factors, rotated_loadings); + } + + show_explained_variance (factor, idata, idata->eval, extracted_eigenvalues, rotated_loadings); factor_matrix_workspace_free (fmw); - show_factor_matrix (factor, idata, factor_matrix); + show_scree (factor, idata); + + show_factor_matrix (factor, idata, + factor->extraction == EXTRACTION_PC ? _("Component Matrix") : _("Factor Matrix"), + factor_matrix); + + if ( factor->rotation != ROT_NONE) + { + show_factor_matrix (factor, idata, + factor->extraction == EXTRACTION_PC ? _("Rotated Component Matrix") : _("Rotated Factor Matrix"), + rotated_factors); + + gsl_matrix_free (rotated_factors); + } + + gsl_vector_free (initial_communalities); gsl_vector_free (extracted_communalities); } + finish: + idata_free (idata); casereader_destroy (r); } + + +