correlation coefficient

correlation coefficient

[Patrick Alken]

Hi,

  Is there any interest in putting a new function in the
statistics area for calculating the Pearson correlation coefficient?
I think this can be done safely in gsl by just using

r = gsl_stats_covariance(x,y) / (gsl_stats_sd(x) * gsl_stats_sd(y))

but it would be more efficient to calculate everything in 1 pass
through the data and I believe there is a stable algorithm to do
this (similar to how the mean/variance is calculated). This is
such a common function for people who work with data so I think
it'd be nice to have it in gsl :)

Patrick Alken

Re: correlation coefficient

[Ben Klemens]

>   Is there any interest in putting a new function in the
> statistics area for calculating the Pearson correlation coefficient?
> I think this can be done safely in gsl by just using
> 
> r = gsl_stats_covariance(x,y) / (gsl_stats_sd(x) * gsl_stats_sd(y))
> 
> but it would be more efficient to calculate everything in 1 pass
> through the data and I believe there is a stable algorithm to do
> this (similar to how the mean/variance is calculated). This is
> such a common function for people who work with data so I think
> it'd be nice to have it in gsl :)

I've been working on a library of stats functions to complement
the GSL, so it naturally includes a correlation matrix
function (apop_correlation_matrix). The library home page is at
http://apophenia.info . There's an accompanying book whose home page
(this week) is at http://avocado.econ.jhu.edu/modeling .

Responding to your request for a covariance with an entire package may
be overkill, but I assume if you're looking for one statistic, you're
probably looking for several more.

Regards,

BK

Re: correlation coefficient

[Brian Gough]

At Thu, 15 Mar 2007 17:06:43 -0600,
Patrick Alken wrote:
>   Is there any interest in putting a new function in the
> statistics area for calculating the Pearson correlation coefficient?
> I think this can be done safely in gsl by just using
> 
> r = gsl_stats_covariance(x,y) / (gsl_stats_sd(x) * gsl_stats_sd(y))
> 
> but it would be more efficient to calculate everything in 1 pass
> through the data and I believe there is a stable algorithm to do
> this (similar to how the mean/variance is calculated). 

Yes, sounds like a good idea to me. Go ahead and add it in
covariance_source.c if you have the 1-pass algorithm.

-- 
Brian Gough

RE: correlation coefficient

[keith.briggs]

>>   Is there any interest in putting a new function in the
>> statistics area for calculating the Pearson correlation coefficient?
>> ... it would be more efficient to calculate everything in 1 pass
>> through the data 
>Yes, sounds like a good idea to me. Go ahead and add it in
>covariance_source.c if you have the 1-pass algorithm.

A one-pass algorithm for this and also for polynomial least-squares can be found at http://keithbriggs.info/pipemath.html.

Keith 

	

Re: correlation coefficient

[James Theiler]

On Fri, 16 Mar 2007, Brian Gough wrote:

] At Thu, 15 Mar 2007 17:06:43 -0600,
] Patrick Alken wrote:
] >   Is there any interest in putting a new function in the
] > statistics area for calculating the Pearson correlation coefficient?
] > I think this can be done safely in gsl by just using
] > 
] > r = gsl_stats_covariance(x,y) / (gsl_stats_sd(x) * gsl_stats_sd(y))
] > 
] > but it would be more efficient to calculate everything in 1 pass
] > through the data and I believe there is a stable algorithm to do
] > this (similar to how the mean/variance is calculated). 
] 
] Yes, sounds like a good idea to me. Go ahead and add it in
] covariance_source.c if you have the 1-pass algorithm.
] 
] 

be sure to include Pearson in the name of the function, since there
are also Spearman's and Kendall's correlation statistics.  (on second
thought, contradicting myself, those two are specialized nonparametric
measures, and so maybe it's reasonable to have Pearson's be the
default.)

jt

-- 
James Theiler
MS-B244, ISR-2, LANL; Los Alamos, NM 87544
Space and Remote Sensing Sciences; Los Alamos National Laboratory
http://public.lanl.gov/jt

Re: correlation coefficient

[Patrick Alken]

> be sure to include Pearson in the name of the function, since there
> are also Spearman's and Kendall's correlation statistics.  (on second
> thought, contradicting myself, those two are specialized nonparametric
> measures, and so maybe it's reasonable to have Pearson's be the
> default.)

matlab and octave both use the name "corrcoef". I don't know about
matlab, but octave has separate functions "spearman" and "kendall"
for the other statistics, so I was thinking that
gsl_stats_correlation would be ok, and later if people want the
other functions they can name them gsl_stats_spearman etc. I'll
make sure the docs clearly state its the Pearson coefficient
for the function.

Patrick Alken

Re: correlation coefficient

[Patrick Alken]

> Yes, sounds like a good idea to me. Go ahead and add it in
> covariance_source.c if you have the 1-pass algorithm.

Ok, added under gsl_stats_correlation. I tested it successfully with
the data set:

x = 9.0e9 + i + 1
y = 9.0e9 - i - 1

for i = 1..100

which normally causes bad codes to fail. Also tested it with
lots of random x/y vectors and compared against the result of:

gsl_stats_covariance(x,y) / (gsl_stats_sd(x) * gsl_stats_sd(y))

all with excellent results - error was below GSL_DBL_EPSILON in all
cases - the test codes in statistics/ aren't really set up for
this type of exhaustive testing...so I just adapted the test that
is in there for the correlation routine.

Patrick Alken

Re: correlation coefficient

[Brian Gough]

At Fri, 16 Mar 2007 16:15:32 -0600,
Patrick Alken wrote:
> Ok, added under gsl_stats_correlation. I tested it successfully with
> the data set...

Looks good.  Incidentally if you didn't find it there is a script
scripts/mkheaders.pl which can be used to update the various headers
automatically from the 'float' version.

-- 
Brian Gough

Re: correlation coefficient

[Patrick Alken]

> Looks good.  Incidentally if you didn't find it there is a script
> scripts/mkheaders.pl which can be used to update the various headers
> automatically from the 'float' version.

Ah..good to know :-)