1 /* PSPP - a program for statistical analysis.
2 Copyright (C) 2009 Free Software Foundation, Inc.
4 This program is free software: you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation, either version 3 of the License, or
7 (at your option) any later version.
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
14 You should have received a copy of the GNU General Public License
15 along with this program. If not, see <http://www.gnu.org/licenses/>. */
19 #include <libpspp/assertion.h>
20 #include "covariance.h"
21 #include <gl/xalloc.h>
23 #include <gsl/gsl_matrix.h>
24 #include <data/case.h>
25 #include <data/variable.h>
26 #include <libpspp/misc.h>
27 #include "categoricals.h"
29 #define n_MOMENTS (MOMENT_VARIANCE + 1)
32 /* Create a new matrix of NEW_SIZE x NEW_SIZE and copy the elements of
33 matrix IN into it. IN must be a square matrix, and in normal usage
34 it will be smaller than NEW_SIZE.
35 IN is destroyed by this function. The return value must be destroyed
36 when no longer required.
39 resize_matrix (gsl_matrix *in, size_t new_size)
43 gsl_matrix *out = NULL;
45 assert (in->size1 == in->size2);
47 if (new_size <= in->size1)
50 out = gsl_matrix_calloc (new_size, new_size);
52 for (i = 0; i < in->size1; ++i)
54 for (j = 0; j < in->size2; ++j)
56 double x = gsl_matrix_get (in, i, j);
58 gsl_matrix_set (out, i, j, x);
69 /* The variables for which the covariance matrix is to be calculated. */
71 const struct variable **vars;
73 /* Categorical variables. */
74 struct categoricals *categoricals;
76 /* Array containing number of categories per categorical variable. */
79 /* Dimension of the covariance matrix. */
82 /* The weight variable (or NULL if none) */
83 const struct variable *wv;
85 /* A set of matrices containing the 0th, 1st and 2nd moments */
88 /* The class of missing values to exclude */
89 enum mv_class exclude;
91 /* An array of doubles representing the covariance matrix.
92 Only the top triangle is included, and no diagonals */
96 /* 1 for single pass algorithm;
97 2 for double pass algorithm
102 0 : No pass has been made
103 1 : First pass has been started
104 2 : Second pass has been
106 IE: How many passes have been (partially) made. */
109 /* Flags indicating that the first case has been seen */
110 bool pass_one_first_case_seen;
111 bool pass_two_first_case_seen;
116 /* Return a matrix containing the M th moments.
117 The matrix is of size NxN where N is the number of variables.
118 Each row represents the moments of a variable.
119 In the absence of missing values, the columns of this matrix will
120 be identical. If missing values are involved, then element (i,j)
121 is the moment of the i th variable, when paired with the j th variable.
124 covariance_moments (const struct covariance *cov, int m)
126 return cov->moments[m];
131 /* Create a covariance struct.
134 covariance_1pass_create (size_t n_vars, const struct variable **vars,
135 const struct variable *weight, enum mv_class exclude)
138 struct covariance *cov = xzalloc (sizeof *cov);
142 cov->pass_one_first_case_seen = cov->pass_two_first_case_seen = false;
147 cov->n_vars = n_vars;
150 cov->moments = xmalloc (sizeof *cov->moments * n_MOMENTS);
152 for (i = 0; i < n_MOMENTS; ++i)
153 cov->moments[i] = gsl_matrix_calloc (n_vars, n_vars);
155 cov->exclude = exclude;
157 cov->n_cm = (n_vars * (n_vars - 1) ) / 2;
160 cov->cm = xcalloc (sizeof *cov->cm, cov->n_cm);
161 cov->categoricals = NULL;
167 Create a covariance struct for a two-pass algorithm. If categorical
168 variables are involed, the dimension cannot be know until after the
169 first data pass, so the actual covariances will not be allocated
173 covariance_2pass_create (size_t n_vars, const struct variable **vars,
174 size_t n_catvars, const struct variable **catvars,
175 const struct variable *wv, enum mv_class exclude)
178 struct covariance *cov = xmalloc (sizeof *cov);
182 cov->pass_one_first_case_seen = cov->pass_two_first_case_seen = false;
187 cov->n_vars = n_vars;
190 cov->moments = xmalloc (sizeof *cov->moments * n_MOMENTS);
192 for (i = 0; i < n_MOMENTS; ++i)
193 cov->moments[i] = gsl_matrix_calloc (n_vars, n_vars);
195 cov->exclude = exclude;
200 cov->categoricals = categoricals_create (catvars, n_catvars, wv, exclude);
205 /* Return an integer, which can be used to index
206 into COV->cm, to obtain the I, J th element
207 of the covariance matrix. If COV->cm does not
208 contain that element, then a negative value
212 cm_idx (const struct covariance *cov, int i, int j)
215 const int n2j = cov->dim - 2 - j;
216 const int nj = cov->dim - 2 ;
219 assert (j < cov->dim);
224 if (j >= cov->dim - 1)
231 as -= n2j * (n2j + 1) ;
239 Returns true iff the variable corresponding to the Ith element of the covariance matrix
240 has a missing value for case C
243 is_missing (const struct covariance *cov, int i, const struct ccase *c)
245 const struct variable *var = i < cov->n_vars ?
247 categoricals_get_variable_by_subscript (cov->categoricals, i - cov->n_vars);
249 const union value *val = case_data (c, var);
251 return var_is_value_missing (var, val, cov->exclude);
256 get_val (const struct covariance *cov, int i, const struct ccase *c)
258 if ( i < cov->n_vars)
260 const struct variable *var = cov->vars[i];
262 const union value *val = case_data (c, var);
267 return categoricals_get_binary_by_subscript (cov->categoricals, i - cov->n_vars, c);
271 dump_matrix (const gsl_matrix *m)
275 for (i = 0 ; i < m->size1; ++i)
277 for (j = 0 ; j < m->size2; ++j)
278 printf ("%02f ", gsl_matrix_get (m, i, j));
283 /* Call this function for every case in the data set */
285 covariance_accumulate_pass1 (struct covariance *cov, const struct ccase *c)
288 const double weight = cov->wv ? case_data (c, cov->wv)->f : 1.0;
290 assert (cov->passes == 2);
291 if (!cov->pass_one_first_case_seen)
293 assert (cov->state == 0);
297 categoricals_update (cov->categoricals, c);
299 for (i = 0 ; i < cov->dim; ++i)
301 double v1 = get_val (cov, i, c);
303 if ( is_missing (cov, i, c))
306 for (j = 0 ; j < cov->dim; ++j)
310 if ( is_missing (cov, j, c))
313 for (m = 0 ; m <= MOMENT_MEAN; ++m)
315 double *x = gsl_matrix_ptr (cov->moments[m], i, j);
323 cov->pass_one_first_case_seen = true;
327 /* Call this function for every case in the data set */
329 covariance_accumulate_pass2 (struct covariance *cov, const struct ccase *c)
332 const double weight = cov->wv ? case_data (c, cov->wv)->f : 1.0;
334 assert (cov->passes == 2);
335 assert (cov->state >= 1);
337 if (! cov->pass_two_first_case_seen)
340 assert (cov->state == 1);
343 cov->dim = cov->n_vars +
344 categoricals_total (cov->categoricals) - categoricals_get_n_variables (cov->categoricals);
346 cov->n_cm = (cov->dim * (cov->dim - 1) ) / 2;
347 cov->cm = xcalloc (sizeof *cov->cm, cov->n_cm);
349 /* Grow the moment matrices so that they're large enough to accommodate the
350 categorical elements */
351 for (i = 0; i < n_MOMENTS; ++i)
353 cov->moments[i] = resize_matrix (cov->moments[i], cov->dim);
356 categoricals_done (cov->categoricals);
358 /* Populate the moments matrices with the categorical value elements */
359 for (i = cov->n_vars; i < cov->dim; ++i)
361 for (j = 0 ; j < cov->dim ; ++j) /* FIXME: This is WRONG !!! */
363 double w = categoricals_get_weight_by_subscript (cov->categoricals, i - cov->n_vars);
365 gsl_matrix_set (cov->moments[MOMENT_NONE], i, j, w);
367 w = categoricals_get_sum_by_subscript (cov->categoricals, i - cov->n_vars);
369 gsl_matrix_set (cov->moments[MOMENT_MEAN], i, j, w);
373 /* FIXME: This is WRONG!! It must be fixed to properly handle missing values. For
374 now it assumes there are none */
375 for (m = 0 ; m < n_MOMENTS; ++m)
377 for (i = 0 ; i < cov->dim ; ++i)
379 double x = gsl_matrix_get (cov->moments[m], i, cov->n_vars -1);
380 for (j = cov->n_vars; j < cov->dim; ++j)
382 gsl_matrix_set (cov->moments[m], i, j, x);
387 /* Divide the means by the number of samples */
388 for (i = 0; i < cov->dim; ++i)
390 for (j = 0; j < cov->dim; ++j)
392 double *x = gsl_matrix_ptr (cov->moments[MOMENT_MEAN], i, j);
393 *x /= gsl_matrix_get (cov->moments[MOMENT_NONE], i, j);
398 for (i = 0 ; i < cov->dim; ++i)
400 double v1 = get_val (cov, i, c);
402 if ( is_missing (cov, i, c))
405 for (j = 0 ; j < cov->dim; ++j)
409 double v2 = get_val (cov, j, c);
411 const double s = pow2 (v1 - gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j)) * weight;
413 if ( is_missing (cov, j, c))
417 double *x = gsl_matrix_ptr (cov->moments[MOMENT_VARIANCE], i, j);
422 (v1 - gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j))
424 (v2 - gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j))
428 idx = cm_idx (cov, i, j);
436 cov->pass_two_first_case_seen = true;
440 /* Call this function for every case in the data set.
441 After all cases have been passed, call covariance_calculate
444 covariance_accumulate (struct covariance *cov, const struct ccase *c)
447 const double weight = cov->wv ? case_data (c, cov->wv)->f : 1.0;
449 assert (cov->passes == 1);
451 if ( !cov->pass_one_first_case_seen)
453 assert ( cov->state == 0);
457 for (i = 0 ; i < cov->dim; ++i)
459 const union value *val1 = case_data (c, cov->vars[i]);
461 if ( is_missing (cov, i, c))
464 for (j = 0 ; j < cov->dim; ++j)
468 const union value *val2 = case_data (c, cov->vars[j]);
470 if ( is_missing (cov, j, c))
473 idx = cm_idx (cov, i, j);
476 cov->cm [idx] += val1->f * val2->f * weight;
479 for (m = 0 ; m < n_MOMENTS; ++m)
481 double *x = gsl_matrix_ptr (cov->moments[m], i, j);
489 cov->pass_one_first_case_seen = true;
494 Allocate and return a gsl_matrix containing the covariances of the
498 cm_to_gsl (struct covariance *cov)
501 gsl_matrix *m = gsl_matrix_calloc (cov->dim, cov->dim);
503 /* Copy the non-diagonal elements from cov->cm */
504 for ( j = 0 ; j < cov->dim - 1; ++j)
506 for (i = j+1 ; i < cov->dim; ++i)
508 double x = cov->cm [cm_idx (cov, i, j)];
509 gsl_matrix_set (m, i, j, x);
510 gsl_matrix_set (m, j, i, x);
514 /* Copy the diagonal elements from cov->moments[2] */
515 for (j = 0 ; j < cov->dim ; ++j)
517 double sigma = gsl_matrix_get (cov->moments[2], j, j);
518 gsl_matrix_set (m, j, j, sigma);
525 static const gsl_matrix *
526 covariance_calculate_double_pass (struct covariance *cov)
529 for (i = 0 ; i < cov->dim; ++i)
531 for (j = 0 ; j < cov->dim; ++j)
534 double *x = gsl_matrix_ptr (cov->moments[MOMENT_VARIANCE], i, j);
535 *x /= gsl_matrix_get (cov->moments[MOMENT_NONE], i, j);
537 idx = cm_idx (cov, i, j);
541 *x /= gsl_matrix_get (cov->moments[MOMENT_NONE], i, j);
546 return cm_to_gsl (cov);
549 static const gsl_matrix *
550 covariance_calculate_single_pass (struct covariance *cov)
555 for (m = 0; m < n_MOMENTS; ++m)
557 /* Divide the moments by the number of samples */
560 for (i = 0 ; i < cov->dim; ++i)
562 for (j = 0 ; j < cov->dim; ++j)
564 double *x = gsl_matrix_ptr (cov->moments[m], i, j);
565 *x /= gsl_matrix_get (cov->moments[0], i, j);
567 if ( m == MOMENT_VARIANCE)
568 *x -= pow2 (gsl_matrix_get (cov->moments[1], i, j));
574 /* Centre the moments */
575 for ( j = 0 ; j < cov->dim - 1; ++j)
577 for (i = j + 1 ; i < cov->dim; ++i)
579 double *x = &cov->cm [cm_idx (cov, i, j)];
581 *x /= gsl_matrix_get (cov->moments[0], i, j);
584 gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j)
586 gsl_matrix_get (cov->moments[MOMENT_MEAN], j, i);
590 return cm_to_gsl (cov);
595 Return a pointer to gsl_matrix containing the pairwise covariances.
596 The matrix remains owned by the COV object, and must not be freed.
597 Call this function only after all data have been accumulated.
600 covariance_calculate (struct covariance *cov)
602 assert ( cov->state > 0 );
607 return covariance_calculate_single_pass (cov);
610 return covariance_calculate_double_pass (cov);
618 Covariance computed without dividing by the sample size.
620 static const gsl_matrix *
621 covariance_calculate_double_pass_unnormalized (struct covariance *cov)
624 for (i = 0 ; i < cov->dim; ++i)
626 for (j = 0 ; j < cov->dim; ++j)
629 double *x = gsl_matrix_ptr (cov->moments[MOMENT_VARIANCE], i, j);
631 idx = cm_idx (cov, i, j);
639 return cm_to_gsl (cov);
642 static const gsl_matrix *
643 covariance_calculate_single_pass_unnormalized (struct covariance *cov)
648 for (i = 0 ; i < cov->dim; ++i)
650 for (j = 0 ; j < cov->dim; ++j)
652 double *x = gsl_matrix_ptr (cov->moments[MOMENT_VARIANCE], i, j);
653 *x -= pow2 (gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j))
654 / gsl_matrix_get (cov->moments[MOMENT_NONE], i, j);
657 for ( j = 0 ; j < cov->dim - 1; ++j)
659 for (i = j + 1 ; i < cov->dim; ++i)
661 double *x = &cov->cm [cm_idx (cov, i, j)];
664 gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j)
666 gsl_matrix_get (cov->moments[MOMENT_MEAN], j, i)
667 / gsl_matrix_get (cov->moments[MOMENT_NONE], i, j);
671 return cm_to_gsl (cov);
676 Return a pointer to gsl_matrix containing the pairwise covariances.
677 The matrix remains owned by the COV object, and must not be freed.
678 Call this function only after all data have been accumulated.
681 covariance_calculate_unnormalized (struct covariance *cov)
683 assert ( cov->state > 0 );
688 return covariance_calculate_single_pass_unnormalized (cov);
691 return covariance_calculate_double_pass_unnormalized (cov);
700 /* Destroy the COV object */
702 covariance_destroy (struct covariance *cov)
706 categoricals_destroy (cov->categoricals);
708 for (i = 0; i < n_MOMENTS; ++i)
709 gsl_matrix_free (cov->moments[i]);