double econverge;
int iterations;
+ double rconverge;
+
/* Format */
double blank;
bool sort;
}
+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);
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;
double d = 0.0;
double x, y;
double phi;
+
for (p = 0; p < normalised->size1; ++p)
{
double jv = gsl_matrix_get (normalised, p, j);
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,
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)
{
*lambda = - *lambda;
}
}
-
- // dump_matrix (result);
}
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_force_match (lexer, ')');
}
}
+ else if (lex_match_id (lexer, "RCONVERGE"))
+ {
+ if ( lex_force_match (lexer, '('))
+ {
+ lex_force_num (lexer);
+ factor.rconverge = lex_number (lexer);
+ lex_get (lexer);
+ lex_force_match (lexer, ')');
+ }
+ }
else if (lex_match_id (lexer, "ITERATE"))
{
if ( lex_force_match (lexer, '('))
}
}
+ if ( factor.rotation == ROT_NONE )
+ factor.print &= ~PRINT_ROTATION;
+
if ( ! run_factor (ds, &factor))
goto error;
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;
double e_total = 0.0;
double e_cum = 0.0;
+ double r_cum = 0.0;
+
int nc = heading_columns;
if (factor->print & PRINT_EXTRACTION)
e_total = i_total;
}
-
for (i = 0 ; i < factor->n_vars; ++i)
{
const double i_lambda = gsl_vector_get (initial_eigenvalues, i);
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)
tab_double (t, c++, i + heading_rows, 0, i_cum, NULL);
}
+
if (factor->print & PRINT_EXTRACTION)
{
if (i < idata->n_extractions)
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);
}
{
+ 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);
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);
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);
}
+
gsl_vector_free (initial_communalities);
gsl_vector_free (extracted_communalities);
}
casereader_destroy (r);
}
+
+
+
projects / pspp-builds.git / blobdiff