X-Git-Url: https://pintos-os.org/cgi-bin/gitweb.cgi?a=blobdiff_plain;f=src%2Flanguage%2Fstats%2Ffactor.c;h=a19cf69e16e9e3f133fc74495fae07f599903e38;hb=81579d9e9f994fb2908f50af41c3eb033d216e58;hp=d3bba77eaa3c68aac1136e4817c958daa8ceefc9;hpb=e88057263e1a2ade0c166d21f5685d6dc2797257;p=pspp-builds.git diff --git a/src/language/stats/factor.c b/src/language/stats/factor.c index d3bba77e..a19cf69e 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, 2010 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,30 +23,25 @@ #include #include -#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/dictionary.h" +#include "data/format.h" +#include "data/procedure.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) @@ -163,6 +157,8 @@ struct cmd_factor double econverge; int iterations; + double rconverge; + /* Format */ double blank; bool sort; @@ -173,7 +169,7 @@ struct idata /* Intermediate values used in calculation */ const gsl_matrix *corr ; /* The correlation matrix */ - const gsl_matrix *cov ; /* The covariance matrix */ + gsl_matrix *cov ; /* The covariance matrix */ const gsl_matrix *n ; /* Matrix of number of samples */ gsl_vector *eval ; /* The eigenvalues */ @@ -204,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) { @@ -247,6 +246,7 @@ dump_vector (const gsl_vector *v) } printf ("\n"); } +#endif static int @@ -539,11 +539,41 @@ clone_matrix (const gsl_matrix *m) } +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 gsl_matrix *unrot, const gsl_vector *communalities, enum rotation_type rot_type, gsl_matrix *result) +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); @@ -571,11 +601,19 @@ rotate (const gsl_matrix *unrot, const gsl_vector *communalities, enum rotation_ gsl_matrix_free (h_sqrt_inv); + /* Now perform the rotation iterations */ - for (i = 0 ; i < 25 ; ++i) + + 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; @@ -585,6 +623,7 @@ rotate (const gsl_matrix *unrot, const gsl_vector *communalities, enum rotation_ double d = 0.0; double x, y; double phi; + for (p = 0; p < normalised->size1; ++p) { double jv = gsl_matrix_get (normalised, p, j); @@ -598,18 +637,38 @@ rotate (const gsl_matrix *unrot, const gsl_vector *communalities, enum rotation_ d += 2 * u * v; } - rotation_coeff [rot_type] (&x, &y, a, b, c, d, normalised); + 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, @@ -618,15 +677,20 @@ rotate (const gsl_matrix *unrot, const gsl_vector *communalities, enum rotation_ gsl_matrix_free (h_sqrt); - /* reflect negative sums */ + /* 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) { @@ -634,8 +698,6 @@ rotate (const gsl_matrix *unrot, const gsl_vector *communalities, enum rotation_ *lambda = - *lambda; } } - - // dump_matrix (result); } @@ -712,21 +774,24 @@ 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)) @@ -735,14 +800,14 @@ cmd_factor (struct lexer *lexer, struct dataset *ds) if (factor.n_vars < 2) msg (MW, _("Factor analysis on a single variable is not useful.")); - while (lex_token (lexer) != '.') + while (lex_token (lexer) != T_ENDCMD) { - lex_match (lexer, '/'); + lex_match (lexer, T_SLASH); 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")) { @@ -762,8 +827,8 @@ cmd_factor (struct lexer *lexer, struct dataset *ds) } 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")) { @@ -782,8 +847,8 @@ cmd_factor (struct lexer *lexer, struct dataset *ds) } 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) { /* VARIMAX and DEFAULT are defaults */ if (lex_match_id (lexer, "VARIMAX") || lex_match_id (lexer, "DEFAULT")) @@ -811,47 +876,57 @@ cmd_factor (struct lexer *lexer, struct dataset *ds) } 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")) @@ -870,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")) { @@ -898,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")) { @@ -907,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")) @@ -930,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")) { @@ -1004,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")) { @@ -1041,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; @@ -1293,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; @@ -1309,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) @@ -1396,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); @@ -1405,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) @@ -1420,6 +1504,7 @@ 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 < idata->n_extractions) @@ -1430,6 +1515,17 @@ 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); @@ -1531,7 +1627,7 @@ show_correlation_matrix (const struct cmd_factor *factor, const struct idata *id 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")); + tab_text (t, 0, y, TAT_TITLE, _("Sig. (1-tailed)")); for (i = 0; i < factor->n_vars; ++i) { @@ -1584,7 +1680,7 @@ do_factor (const struct cmd_factor *factor, struct casereader *r) const gsl_matrix *analysis_matrix; 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)) @@ -1594,6 +1690,12 @@ do_factor (const struct cmd_factor *factor, struct casereader *r) 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); idata->n = covariance_moments (cov, MOMENT_NONE); @@ -1686,6 +1788,9 @@ do_factor (const struct cmd_factor *factor, struct casereader *r) } { + 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); @@ -1745,7 +1850,16 @@ do_factor (const struct cmd_factor *factor, struct casereader *r) 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); @@ -1757,10 +1871,6 @@ do_factor (const struct cmd_factor *factor, struct casereader *r) if ( factor->rotation != ROT_NONE) { - gsl_matrix *rotated_factors = gsl_matrix_calloc (factor_matrix->size1, factor_matrix->size2); - - rotate (factor_matrix, extracted_communalities, factor->rotation, rotated_factors); - show_factor_matrix (factor, idata, factor->extraction == EXTRACTION_PC ? _("Rotated Component Matrix") : _("Rotated Factor Matrix"), rotated_factors); @@ -1769,6 +1879,7 @@ do_factor (const struct cmd_factor *factor, struct casereader *r) } + gsl_vector_free (initial_communalities); gsl_vector_free (extracted_communalities); } @@ -1779,3 +1890,6 @@ do_factor (const struct cmd_factor *factor, struct casereader *r) casereader_destroy (r); } + + +