/gsl/trunk

To get this branch, use:
bzr branch /lh/gsl/trunk
2194 by bjg
added covariance function
1
/* statistics/covar_source.c
2
 * 
4276 by bjg
updated copyright years
3
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Jim Davies, Brian Gough
2194 by bjg
added covariance function
4
 * 
5
 * This program is free software; you can redistribute it and/or modify
6
 * it under the terms of the GNU General Public License as published by
4212 by bjg
updated to gpl v3
7
 * the Free Software Foundation; either version 3 of the License, or (at
2194 by bjg
added covariance function
8
 * your option) any later version.
9
 * 
10
 * This program is distributed in the hope that it will be useful, but
11
 * WITHOUT ANY WARRANTY; without even the implied warranty of
12
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
13
 * General Public License for more details.
14
 * 
15
 * You should have received a copy of the GNU General Public License
16
 * along with this program; if not, write to the Free Software
3825 by bjg
updated FSF address in copyright notices
17
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
2194 by bjg
added covariance function
18
 */
19
20
static double
21
FUNCTION(compute,covariance) (const BASE data1[], const size_t stride1,
22
                              const BASE data2[], const size_t stride2,
23
                              const size_t n, 
24
                              const double mean1, const double mean2);
25
26
static double
27
FUNCTION(compute,covariance) (const BASE data1[], const size_t stride1,
28
                              const BASE data2[], const size_t stride2,
29
                              const size_t n, 
30
                              const double mean1, const double mean2)
31
{
32
  /* takes a dataset and finds the covariance */
33
34
  long double covariance = 0 ;
35
36
  size_t i;
37
38
  /* find the sum of the squares */
39
  for (i = 0; i < n; i++)
40
    {
41
      const long double delta1 = (data1[i * stride1] - mean1);
42
      const long double delta2 = (data2[i * stride2] - mean2);
43
      covariance += (delta1 * delta2 - covariance) / (i + 1);
44
    }
45
46
  return covariance ;
47
}
48
49
double 
50
FUNCTION(gsl_stats,covariance_m) (const BASE data1[], const size_t stride1, 
51
                                  const BASE data2[], const size_t stride2, 
52
                                  const size_t n, 
53
                                  const double mean1, const double mean2)
54
{
55
  const double covariance = FUNCTION(compute,covariance) (data1, stride1,
56
                                                          data2, stride2,
57
                                                          n, 
58
                                                          mean1, mean2);
59
  
60
  return covariance * ((double)n / (double)(n - 1));
61
}
62
63
double 
64
FUNCTION(gsl_stats,covariance) (const BASE data1[], const size_t stride1,
65
                                const BASE data2[], const size_t stride2,
66
                                const size_t n)
67
{
68
  const double mean1 = FUNCTION(gsl_stats,mean) (data1, stride1, n);
69
  const double mean2 = FUNCTION(gsl_stats,mean) (data2, stride2, n);
70
71
  return FUNCTION(gsl_stats,covariance_m)(data1, stride1, 
72
                                          data2, stride2, 
73
                                          n, 
74
                                          mean1, mean2);
75
}
76
4152 by palken
added function gsl_stats_correlation
77
/*
78
gsl_stats_correlation()
79
  Calculate Pearson correlation = cov(X, Y) / (sigma_X * sigma_Y)
80
This routine efficiently computes the correlation in one pass of the
81
data and makes use of the algorithm described in:
82
83
B. P. Welford, "Note on a Method for Calculating Corrected Sums of
84
Squares and Products", Technometrics, Vol 4, No 3, 1962.
85
86
This paper derives a numerically stable recurrence to compute a sum
87
of products
88
89
S = sum_{i=1..N} [ (x_i - mu_x) * (y_i - mu_y) ]
90
91
with the relation
92
93
S_n = S_{n-1} + ((n-1)/n) * (x_n - mu_x_{n-1}) * (y_n - mu_y_{n-1})
94
*/
95
96
double
97
FUNCTION(gsl_stats,correlation) (const BASE data1[], const size_t stride1,
98
                                 const BASE data2[], const size_t stride2,
99
                                 const size_t n)
100
{
101
  size_t i;
102
  long double sum_xsq = 0.0;
103
  long double sum_ysq = 0.0;
104
  long double sum_cross = 0.0;
105
  long double ratio;
106
  long double delta_x, delta_y;
107
  long double mean_x, mean_y;
108
  long double r;
109
110
  /*
111
   * Compute:
112
   * sum_xsq = Sum [ (x_i - mu_x)^2 ],
113
   * sum_ysq = Sum [ (y_i - mu_y)^2 ] and
114
   * sum_cross = Sum [ (x_i - mu_x) * (y_i - mu_y) ]
115
   * using the above relation from Welford's paper
116
   */
117
118
  mean_x = data1[0 * stride1];
119
  mean_y = data2[0 * stride2];
120
121
  for (i = 1; i < n; ++i)
122
    {
123
      ratio = i / (i + 1.0);
124
      delta_x = data1[i * stride1] - mean_x;
125
      delta_y = data2[i * stride2] - mean_y;
126
      sum_xsq += delta_x * delta_x * ratio;
127
      sum_ysq += delta_y * delta_y * ratio;
128
      sum_cross += delta_x * delta_y * ratio;
129
      mean_x += delta_x / (i + 1.0);
130
      mean_y += delta_y / (i + 1.0);
131
    }
132
133
  r = sum_cross / (sqrt(sum_xsq) * sqrt(sum_ysq));
134
135
  return r;
136
}
4816 by Patrick Alken
added gsl_stats_spearman function
137
4818 by Patrick Alken
cleanup
138
/*
139
gsl_stats_spearman()
140
  Compute Spearman rank correlation coefficient
141
142
Inputs: data1   - data1 vector
143
        stride1 - stride of data1
144
        data2   - data2 vector
145
        stride2 - stride of data2
146
        n       - number of elements in data1 and data2
147
        work    - additional workspace of size 2*n
148
149
Return: Spearman rank correlation coefficient
150
*/
151
4816 by Patrick Alken
added gsl_stats_spearman function
152
double
153
FUNCTION(gsl_stats,spearman) (const BASE data1[], const size_t stride1,
154
                              const BASE data2[], const size_t stride2,
4818 by Patrick Alken
cleanup
155
                              const size_t n, double work[])
4816 by Patrick Alken
added gsl_stats_spearman function
156
{
157
  size_t i;
158
  gsl_vector_view ranks1 = gsl_vector_view_array(&work[0], n);
159
  gsl_vector_view ranks2 = gsl_vector_view_array(&work[n], n);
160
  double r;
161
162
  for (i = 0; i < n; ++i)
163
    {
164
      gsl_vector_set(&ranks1.vector, i, data1[i * stride1]);
165
      gsl_vector_set(&ranks2.vector, i, data2[i * stride2]);
166
    }
167
168
  /* sort data1 and update data2 at same time; compute rank of data1 */
169
  gsl_sort_vector2(&ranks1.vector, &ranks2.vector);
170
  compute_rank(&ranks1.vector);
171
172
  /* now sort data2, updating ranks1 appropriately; compute rank of data2 */
173
  gsl_sort_vector2(&ranks2.vector, &ranks1.vector);
174
  compute_rank(&ranks2.vector);
175
176
  /* compute correlation of rank vectors in double precision */
177
  r = gsl_stats_correlation(ranks1.vector.data, ranks1.vector.stride,
178
                            ranks2.vector.data, ranks2.vector.stride,
179
                            n);
180
181
  return r;
182
}