X-Git-Url: https://pintos-os.org/cgi-bin/gitweb.cgi?a=blobdiff_plain;f=src%2Flanguage%2Fstats%2Ffactor.c;h=fa9987332fb93bdbdfad55a6f3b58e2750e114a6;hb=e2da62d735c597afeef2e0e9b36e5a4a83d7da94;hp=e09eaae8b01a6d8d00df277791ccd4adcf2a0eba;hpb=691c25e36fd1ee722dd35419d6110e3876b99f9c;p=pspp diff --git a/src/language/stats/factor.c b/src/language/stats/factor.c index e09eaae8b0..fa9987332f 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, 2012 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,38 +16,34 @@ #include - #include #include #include #include #include #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/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/cast.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) @@ -156,12 +152,13 @@ struct cmd_factor enum extraction_method extraction; enum plot_opts plot; enum rotation_type rotation; + int rotation_iterations; /* Extraction Criteria */ int n_factors; double min_eigen; double econverge; - int iterations; + int extraction_iterations; double rconverge; @@ -175,7 +172,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 */ @@ -184,6 +181,8 @@ struct idata int n_extractions; gsl_vector *msr ; /* Multiple Squared Regressions */ + + double detR; /* The determinant of the correlation matrix */ }; static struct idata * @@ -206,11 +205,63 @@ 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); + if (id->corr != NULL) + gsl_matrix_free (CONST_CAST (gsl_matrix *, id->corr)); free (id); } +static gsl_matrix * +anti_image (const gsl_matrix *m) +{ + int i, j; + gsl_matrix *a; + assert (m->size1 == m->size2); + + a = gsl_matrix_alloc (m->size1, m->size2); + + for (i = 0; i < m->size1; ++i) + { + for (j = 0; j < m->size2; ++j) + { + double *p = gsl_matrix_ptr (a, i, j); + *p = gsl_matrix_get (m, i, j); + *p /= gsl_matrix_get (m, i, i); + *p /= gsl_matrix_get (m, j, j); + } + } + + return a; +} + + +/* Return the sum of all the elements excluding row N */ +static double +ssq_od_n (const gsl_matrix *m, int n) +{ + int i, j; + double ss = 0; + assert (m->size1 == m->size2); + + assert (n < m->size1); + + for (i = 0; i < m->size1; ++i) + { + if (i == n ) continue; + for (j = 0; j < m->size2; ++j) + { + ss += pow2 (gsl_matrix_get (m, i, j)); + } + } + + return ss; +} + + + #if 0 static void dump_matrix (const gsl_matrix *m) @@ -225,7 +276,6 @@ dump_matrix (const gsl_matrix *m) } } - static void dump_matrix_permute (const gsl_matrix *m, const gsl_permutation *p) { @@ -609,7 +659,7 @@ rotate (const struct cmd_factor *cf, const gsl_matrix *unrot, /* Now perform the rotation iterations */ prev_sv = initial_sv (normalised); - for (i = 0 ; i < cf->iterations ; ++i) + for (i = 0 ; i < cf->rotation_iterations ; ++i) { double sv = 0.0; for (j = 0 ; j < normalised->size2; ++j) @@ -679,6 +729,7 @@ rotate (const struct cmd_factor *cf, const gsl_matrix *unrot, h_sqrt, normalised, 0.0, result); gsl_matrix_free (h_sqrt); + gsl_matrix_free (normalised); /* reflect negative sums and populate the rotated loadings vector*/ @@ -763,9 +814,8 @@ static bool run_factor (struct dataset *ds, const struct cmd_factor *factor); int cmd_factor (struct lexer *lexer, struct dataset *ds) { - bool extraction_seen = false; const struct dictionary *dict = dataset_dict (ds); - + int n_iterations = 25; struct cmd_factor factor; factor.n_vars = 0; factor.vars = NULL; @@ -776,7 +826,8 @@ cmd_factor (struct lexer *lexer, struct dataset *ds) factor.extraction = EXTRACTION_PC; factor.n_factors = 0; factor.min_eigen = SYSMIS; - factor.iterations = 25; + factor.extraction_iterations = 25; + factor.rotation_iterations = 25; factor.econverge = 0.001; factor.blank = 0; @@ -877,6 +928,7 @@ cmd_factor (struct lexer *lexer, struct dataset *ds) goto error; } } + factor.rotation_iterations = n_iterations; } else if (lex_match_id (lexer, "CRITERIA")) { @@ -928,7 +980,7 @@ cmd_factor (struct lexer *lexer, struct dataset *ds) if ( lex_force_match (lexer, T_LPAREN)) { lex_force_int (lexer); - factor.iterations = lex_integer (lexer); + n_iterations = lex_integer (lexer); lex_get (lexer); lex_force_match (lexer, T_RPAREN); } @@ -937,7 +989,7 @@ cmd_factor (struct lexer *lexer, struct dataset *ds) { factor.n_factors = 0; factor.min_eigen = 1; - factor.iterations = 25; + n_iterations = 25; } else { @@ -948,7 +1000,6 @@ cmd_factor (struct lexer *lexer, struct dataset *ds) } else if (lex_match_id (lexer, "EXTRACTION")) { - extraction_seen = true; lex_match (lexer, T_EQUALS); while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH) { @@ -974,6 +1025,7 @@ cmd_factor (struct lexer *lexer, struct dataset *ds) goto error; } } + factor.extraction_iterations = n_iterations; } else if (lex_match_id (lexer, "FORMAT")) { @@ -1053,10 +1105,11 @@ cmd_factor (struct lexer *lexer, struct dataset *ds) { factor.print |= PRINT_INITIAL; } -#if FACTOR_FULLY_IMPLEMENTED else if (lex_match_id (lexer, "KMO")) { + factor.print |= PRINT_KMO; } +#if FACTOR_FULLY_IMPLEMENTED else if (lex_match_id (lexer, "REPR")) { } @@ -1301,7 +1354,7 @@ show_factor_matrix (const struct cmd_factor *factor, struct idata *idata, const tab_title (t, _("Factor Matrix")); */ - tab_title (t, title); + tab_title (t, "%s", title); tab_headers (t, heading_columns, 0, heading_rows, 0); @@ -1489,16 +1542,12 @@ 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); + tab_text_format (t, c++, i + heading_rows, TAB_LEFT | TAT_TITLE, _("%zu"), i + 1); i_cum += i_percent; e_cum += e_percent; - r_cum += r_percent; /* Initial Eigenvalues */ if (factor->print & PRINT_INITIAL) @@ -1520,16 +1569,22 @@ show_explained_variance (const struct cmd_factor * factor, struct idata *idata, } } - 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); - } - } + if (rotated_loadings != NULL) + { + const double r_lambda = gsl_vector_get (rotated_loadings, i); + double r_percent = 100.0 * r_lambda / e_total ; + if (factor->print & PRINT_ROTATION) + { + if (i < idata->n_extractions) + { + r_cum += r_percent; + 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); @@ -1651,22 +1706,9 @@ show_correlation_matrix (const struct cmd_factor *factor, const struct idata *id 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_double (t, 1, nr, 0, idata->detR, NULL); } tab_submit (t); @@ -1697,27 +1739,46 @@ do_factor (const struct cmd_factor *factor, struct casereader *r) if (idata->cov == NULL) { msg (MW, _("The dataset contains no complete observations. No analysis will be performed.")); + covariance_destroy (cov); goto finish; } var_matrix = covariance_moments (cov, MOMENT_VARIANCE); mean_matrix = covariance_moments (cov, MOMENT_MEAN); idata->n = covariance_moments (cov, MOMENT_NONE); + if ( factor->method == METHOD_CORR) { idata->corr = correlation_from_covariance (idata->cov, var_matrix); + analysis_matrix = idata->corr; } else analysis_matrix = idata->cov; + + if (factor->print & PRINT_DETERMINANT + || factor->print & PRINT_KMO) + { + int sign = 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); + idata->detR = gsl_linalg_LU_det (tmp, sign); + gsl_permutation_free (p); + gsl_matrix_free (tmp); + } + if ( factor->print & PRINT_UNIVARIATE) { + const struct fmt_spec *wfmt = factor->wv ? var_get_print_format (factor->wv) : & F_8_0; const int nc = 4; int i; - const struct fmt_spec *wfmt = factor->wv ? var_get_print_format (factor->wv) : & F_8_0; - const int heading_columns = 1; const int heading_rows = 1; @@ -1763,19 +1824,100 @@ do_factor (const struct cmd_factor *factor, struct casereader *r) tab_submit (t); } + if (factor->print & PRINT_KMO) + { + int i; + double sum_ssq_r = 0; + double sum_ssq_a = 0; + + double df = factor->n_vars * ( factor->n_vars - 1) / 2; + + double w = 0; + + + double xsq; + + const int heading_columns = 2; + const int heading_rows = 0; + + const int nr = heading_rows + 4; + const int nc = heading_columns + 1; + + gsl_matrix *a, *x; + + struct tab_table *t = tab_create (nc, nr); + tab_title (t, _("KMO and Bartlett's Test")); + + x = clone_matrix (idata->corr); + gsl_linalg_cholesky_decomp (x); + gsl_linalg_cholesky_invert (x); + + a = anti_image (x); + + for (i = 0; i < x->size1; ++i) + { + sum_ssq_r += ssq_od_n (x, i); + sum_ssq_a += ssq_od_n (a, i); + } + + gsl_matrix_free (a); + gsl_matrix_free (x); + + tab_headers (t, heading_columns, 0, heading_rows, 0); + + /* Outline the box */ + tab_box (t, + TAL_2, TAL_2, + -1, -1, + 0, 0, + nc - 1, nr - 1); + + tab_vline (t, TAL_2, heading_columns, 0, nr - 1); + + tab_text (t, 0, 0, TAT_TITLE | TAB_LEFT, _("Kaiser-Meyer-Olkin Measure of Sampling Adequacy")); + + tab_double (t, 2, 0, 0, sum_ssq_r / (sum_ssq_r + sum_ssq_a), NULL); + + tab_text (t, 0, 1, TAT_TITLE | TAB_LEFT, _("Bartlett's Test of Sphericity")); + + tab_text (t, 1, 1, TAT_TITLE, _("Approx. Chi-Square")); + tab_text (t, 1, 2, TAT_TITLE, _("df")); + tab_text (t, 1, 3, TAT_TITLE, _("Sig.")); + + + /* The literature doesn't say what to do for the value of W when + missing values are involved. The best thing I can think of + is to take the mean average. */ + w = 0; + for (i = 0; i < idata->n->size1; ++i) + w += gsl_matrix_get (idata->n, i, i); + w /= idata->n->size1; + + xsq = w - 1 - (2 * factor->n_vars + 5) / 6.0; + xsq *= -log (idata->detR); + + tab_double (t, 2, 1, 0, xsq, NULL); + tab_double (t, 2, 2, 0, df, &F_8_0); + tab_double (t, 2, 3, 0, gsl_cdf_chisq_Q (xsq, df), NULL); + + + tab_submit (t); + } + show_correlation_matrix (factor, idata); + covariance_destroy (cov); -#if 1 { + gsl_matrix *am = matrix_dup (analysis_matrix); gsl_eigen_symmv_workspace *workspace = gsl_eigen_symmv_alloc (factor->n_vars); - gsl_eigen_symmv (matrix_dup (analysis_matrix), idata->eval, idata->evec, workspace); + gsl_eigen_symmv (am, idata->eval, idata->evec, workspace); gsl_eigen_symmv_free (workspace); + gsl_matrix_free (am); } gsl_eigen_symmv_sort (idata->eval, idata->evec, GSL_EIGEN_SORT_ABS_DESC); -#endif idata->n_extractions = n_extracted_factors (factor, idata); @@ -1817,7 +1959,7 @@ do_factor (const struct cmd_factor *factor, struct casereader *r) gsl_vector_memcpy (initial_communalities, idata->msr); - for (i = 0; i < factor->iterations; ++i) + for (i = 0; i < factor->extraction_iterations; ++i) { double min, max; gsl_vector_memcpy (diff, idata->msr); @@ -1884,6 +2026,8 @@ do_factor (const struct cmd_factor *factor, struct casereader *r) + gsl_matrix_free (factor_matrix); + gsl_vector_free (rotated_loadings); gsl_vector_free (initial_communalities); gsl_vector_free (extracted_communalities); }