}
/*
- Update the covariance matrix with the new entries, assuming that V1
- is categorical and V2 is numeric.
+ Update the covariance matrix with the new entries, assuming that ROW
+ corresponds to a categorical variable and V2 is numeric.
*/
static void
covariance_update_categorical_numeric (struct design_matrix *cov, double mean,
- double weight, double ssize, const struct variable *v1,
- const struct variable *v2, const union value *val1, const union value *val2)
+ size_t row,
+ const struct variable *v2, double x, const union value *val2)
{
- double x;
- size_t i;
size_t col;
- size_t row;
+ double tmp;
- assert (var_is_alpha (v1));
assert (var_is_numeric (v2));
- row = design_matrix_var_to_column (cov, v1);
col = design_matrix_var_to_column (cov, v2);
- for (i = 0; i < cat_get_n_categories (v1); i++)
+ assert (val2 != NULL);
+ tmp = gsl_matrix_get (cov->m, row, col);
+ gsl_matrix_set (cov->m, row, col, (val2->f - mean) * x + tmp);
+ gsl_matrix_set (cov->m, col, row, (val2->f - mean) * x + tmp);
+}
+static void
+column_iterate (struct design_matrix *cov, const struct variable *v,
+ double ssize, double x, const union value *val1, size_t row)
+{
+ size_t col;
+ size_t i;
+ double y;
+ double tmp;
+ const union value *tmp_val;
+
+ col = design_matrix_var_to_column (cov, v);
+ for (i = 0; i < cat_get_n_categories (v) - 1; i++)
{
- row += i;
- x = -1.0 * cat_get_n_categories (v1) / ssize;
- if (i == cat_value_find (v1, val1))
+ col += i;
+ y = -1.0 * cat_get_category_count (i, v) / ssize;
+ tmp_val = cat_subscript_to_value (i, v);
+ if (compare_values (tmp_val, val1, v))
{
- x += 1.0;
+ y += -1.0;
}
- assert (val2 != NULL);
- gsl_matrix_set (cov->m, row, col, (val2->f - mean) * x * weight);
+ tmp = gsl_matrix_get (cov->m, row, col);
+ gsl_matrix_set (cov->m, row, col, x * y + tmp);
+ gsl_matrix_set (cov->m, col, row, x * y + tmp);
}
}
-
/*
- Call this function in the first data pass. The central moments are
+ Call this function in the second data pass. The central moments are
MEAN1 and MEAN2. Any categorical variables should already have their
values summarized in in its OBS_VALS element.
*/
-void covariance_pass_one (struct design_matrix *cov, double mean1, double mean2,
- double weight, double ssize, const struct variable *v1,
+void covariance_pass_two (struct design_matrix *cov, double mean1, double mean2,
+ double ssize, const struct variable *v1,
const struct variable *v2, const union value *val1, const union value *val2)
{
size_t row;
size_t col;
size_t i;
double x;
- double y;
+ const union value *tmp_val;
if (var_is_alpha (v1))
{
- if (var_is_numeric (v2))
- {
- covariance_update_categorical_numeric (cov, mean2, weight, ssize, v1,
- v2, val1, val2);
- }
- else
+ row = design_matrix_var_to_column (cov, v1);
+ for (i = 0; i < cat_get_n_categories (v1) - 1; i++)
{
- row = design_matrix_var_to_column (cov, v1);
- col = design_matrix_var_to_column (cov, v2);
- for (i = 0; i < cat_get_n_categories (v2); i++)
+ row += i;
+ x = -1.0 * cat_get_category_count (i, v1) / ssize;
+ tmp_val = cat_subscript_to_value (i, v1);
+ if (compare_values (tmp_val, val1, v1))
{
- col += i;
- y = -1.0 * cat_get_n_categories (v2) / ssize;
- if (i == cat_value_find (v2, val2))
- {
- y += 1.0;
- }
- gsl_matrix_set (cov->m, row, col, x * y * weight);
- gsl_matrix_set (cov->m, col, row, x * y * weight);
+ x += 1.0;
+ }
+ if (var_is_numeric (v2))
+ {
+ covariance_update_categorical_numeric (cov, mean2, row,
+ v2, x, val2);
+ }
+ else
+ {
+ column_iterate (cov, v1, ssize, x, val1, row);
+ column_iterate (cov, v2, ssize, x, val2, row);
}
}
}
else if (var_is_alpha (v2))
{
- covariance_update_categorical_numeric (cov, mean1, weight, ssize, v2,
- v1, val2, val1);
+ /*
+ Reverse the orders of V1, V2, etc. and put ourselves back
+ in the previous IF scope.
+ */
+ covariance_pass_two (cov, mean2, mean1, ssize, v2, v1, val2, val1);
}
else
{
*/
row = design_matrix_var_to_column (cov, v1);
col = design_matrix_var_to_column (cov, v2);
- x = (val1->f - mean1) * (val2->f - mean2) * weight;
+ x = (val1->f - mean1) * (val2->f - mean2);
+ x += gsl_matrix_get (cov->m, col, row);
gsl_matrix_set (cov->m, row, col, x);
gsl_matrix_set (cov->m, col, row, x);
}