gsl_stats_wcovariance

gsl_stats_wcovariance

[James Bergstra]

[-- Attachment #1: Type: text/plain, Size: 1189 bytes --]


This email concerns the functions

- gsl_stats_wcovariance
- gsl_stats_wcovariance_w

The file wcovar_source.c was in the cvs tree, but was not included in the build
and it appeared to be incomplete.  Under the assumption that it had simply been
forgotten (copyright notice was 2000!), I wrote it up.

I'd appreciate it if someone who knows about unbiased estimators of weighted
sample statistics would have a glance at my algorithm.  I haven't a clue-- I
implemented a straightforward extension of the wvariance routines.

Please find attached (if the mailing list supports attachments) the files that I
modified.  They are from the statistics subdirectory of today's cvs.  In
addition to rewriting wcovar_source.c, I added wcovariance declarations to
headers (double, float, long double) and added wcovar.c to Makefile.am.  Also
attached is a very simple test program.  This program's output agrees with the
cov.wt routine from R.

If anyone encourages me, I'll try to figure out where/how to add the necessary
documentation.  But before I go further, I'd be curious to know why the original
implementation was abandoned.

-- 
James Bergstra
http://www-etud.iro.umontreal.ca/~bergstrj


[-- Attachment #2: gsl_statistics_double.h --]
[-- Type: text/plain, Size: 6457 bytes --]

/* statistics/gsl_statistics_double.h
 * 
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Jim Davies, Brian Gough
 * 
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#ifndef __GSL_STATISTICS_DOUBLE_H__
#define __GSL_STATISTICS_DOUBLE_H__

#include <stddef.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

double gsl_stats_mean (const double data[], const size_t stride, const size_t n);
double gsl_stats_variance (const double data[], const size_t stride, const size_t n);
double gsl_stats_sd (const double data[], const size_t stride, const size_t n);
double gsl_stats_variance_with_fixed_mean (const double data[], const size_t stride, const size_t n, const double mean);
double gsl_stats_sd_with_fixed_mean (const double data[], const size_t stride, const size_t n, const double mean);
double gsl_stats_absdev (const double data[], const size_t stride, const size_t n);
double gsl_stats_skew (const double data[], const size_t stride, const size_t n);
double gsl_stats_kurtosis (const double data[], const size_t stride, const size_t n);
double gsl_stats_lag1_autocorrelation (const double data[], const size_t stride, const size_t n);

double gsl_stats_covariance (const double data1[], const size_t stride1,const double data2[], const size_t stride2, const size_t n);

double gsl_stats_variance_m (const double data[], const size_t stride, const size_t n, const double mean);
double gsl_stats_sd_m (const double data[], const size_t stride, const size_t n, const double mean);
double gsl_stats_absdev_m (const double data[], const size_t stride, const size_t n, const double mean);
double gsl_stats_skew_m_sd (const double data[], const size_t stride, const size_t n, const double mean, const double sd);
double gsl_stats_kurtosis_m_sd (const double data[], const size_t stride, const size_t n, const double mean, const double sd);
double gsl_stats_lag1_autocorrelation_m (const double data[], const size_t stride, const size_t n, const double mean);

double gsl_stats_covariance_m (const double data1[], const size_t stride1,const double data2[], const size_t stride2, const size_t n, const double mean1, const double mean2);

/* DEFINED FOR FLOATING POINT TYPES ONLY */

double gsl_stats_wmean (const double w[], const size_t wstride, const double data[], const size_t stride, const size_t n);
double gsl_stats_wvariance (const double w[], const size_t wstride, const double data[], const size_t stride, const size_t n);
double gsl_stats_wsd (const double w[], const size_t wstride, const double data[], const size_t stride, const size_t n);
double gsl_stats_wvariance_with_fixed_mean (const double w[], const size_t wstride, const double data[], const size_t stride, const size_t n, const double mean);
double gsl_stats_wsd_with_fixed_mean (const double w[], const size_t wstride, const double data[], const size_t stride, const size_t n, const double mean);
double gsl_stats_wabsdev (const double w[], const size_t wstride, const double data[], const size_t stride, const size_t n);
double gsl_stats_wskew (const double w[], const size_t wstride, const double data[], const size_t stride, const size_t n);
double gsl_stats_wkurtosis (const double w[], const size_t wstride, const double data[], const size_t stride, const size_t n);

double gsl_stats_wcovariance (const double w[], const size_t wstride, const double data1[], const size_t stride1,const double data2[], const size_t stride2, const size_t n);

double gsl_stats_wvariance_m (const double w[], const size_t wstride, const double data[], const size_t stride, const size_t n, const double wmean);
double gsl_stats_wsd_m (const double w[], const size_t wstride, const double data[], const size_t stride, const size_t n, const double wmean);
double gsl_stats_wabsdev_m (const double w[], const size_t wstride, const double data[], const size_t stride, const size_t n, const double wmean);
double gsl_stats_wskew_m_sd (const double w[], const size_t wstride, const double data[], const size_t stride, const size_t n, const double wmean, const double wsd);
double gsl_stats_wkurtosis_m_sd (const double w[], const size_t wstride, const double data[], const size_t stride, const size_t n, const double wmean, const double wsd);

double gsl_stats_wcovariance_m (const double w[], const size_t wstride, const double data1[], const size_t stride1,const double data2[], const size_t stride2, const size_t n, const double mean1, const double mean2);

/* END OF FLOATING POINT TYPES */

double gsl_stats_pvariance (const double data1[], const size_t stride1, const size_t n1, const double data2[], const size_t stride2, const size_t n2);
double gsl_stats_ttest (const double data1[], const size_t stride1, const size_t n1, const double data2[], const size_t stride2, const size_t n2);

double gsl_stats_max (const double data[], const size_t stride, const size_t n);
double gsl_stats_min (const double data[], const size_t stride, const size_t n);
void gsl_stats_minmax (double * min, double * max, const double data[], const size_t stride, const size_t n);

size_t gsl_stats_max_index (const double data[], const size_t stride, const size_t n);
size_t gsl_stats_min_index (const double data[], const size_t stride, const size_t n);
void gsl_stats_minmax_index (size_t * min_index, size_t * max_index, const double data[], const size_t stride, const size_t n);

double gsl_stats_median_from_sorted_data (const double sorted_data[], const size_t stride, const size_t n) ;
double gsl_stats_quantile_from_sorted_data (const double sorted_data[], const size_t stride, const size_t n, const double f) ;

__END_DECLS

#endif /* __GSL_STATISTICS_DOUBLE_H__ */

[-- Attachment #3: gsl_statistics_long_double.h --]
[-- Type: text/plain, Size: 7316 bytes --]

/* statistics/gsl_statistics_long_double.h
 * 
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Jim Davies, Brian Gough
 * 
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#ifndef __GSL_STATISTICS_LONG_DOUBLE_H__
#define __GSL_STATISTICS_LONG_DOUBLE_H__

#include <stddef.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

double gsl_stats_long_double_mean (const long double data[], const size_t stride, const size_t n);
double gsl_stats_long_double_variance (const long double data[], const size_t stride, const size_t n);
double gsl_stats_long_double_sd (const long double data[], const size_t stride, const size_t n);
double gsl_stats_long_double_variance_with_fixed_mean (const long double data[], const size_t stride, const size_t n, const double mean);
double gsl_stats_long_double_sd_with_fixed_mean (const long double data[], const size_t stride, const size_t n, const double mean);
double gsl_stats_long_double_absdev (const long double data[], const size_t stride, const size_t n);
double gsl_stats_long_double_skew (const long double data[], const size_t stride, const size_t n);
double gsl_stats_long_double_kurtosis (const long double data[], const size_t stride, const size_t n);
double gsl_stats_long_double_lag1_autocorrelation (const long double data[], const size_t stride, const size_t n);

double gsl_stats_long_double_covariance (const long double data1[], const size_t stride1,const long double data2[], const size_t stride2, const size_t n);

double gsl_stats_long_double_variance_m (const long double data[], const size_t stride, const size_t n, const double mean);
double gsl_stats_long_double_sd_m (const long double data[], const size_t stride, const size_t n, const double mean);
double gsl_stats_long_double_absdev_m (const long double data[], const size_t stride, const size_t n, const double mean);
double gsl_stats_long_double_skew_m_sd (const long double data[], const size_t stride, const size_t n, const double mean, const double sd);
double gsl_stats_long_double_kurtosis_m_sd (const long double data[], const size_t stride, const size_t n, const double mean, const double sd);
double gsl_stats_long_double_lag1_autocorrelation_m (const long double data[], const size_t stride, const size_t n, const double mean);

double gsl_stats_long_double_covariance_m (const long double data1[], const size_t stride1,const long double data2[], const size_t stride2, const size_t n, const double mean1, const double mean2);

/* DEFINED FOR FLOATING POINT TYPES ONLY */

double gsl_stats_long_double_wmean (const long double w[], const size_t wstride, const long double data[], const size_t stride, const size_t n);
double gsl_stats_long_double_wvariance (const long double w[], const size_t wstride, const long double data[], const size_t stride, const size_t n);
double gsl_stats_long_double_wsd (const long double w[], const size_t wstride, const long double data[], const size_t stride, const size_t n);
double gsl_stats_long_double_wvariance_with_fixed_mean (const long double w[], const size_t wstride, const long double data[], const size_t stride, const size_t n, const double mean);
double gsl_stats_long_double_wsd_with_fixed_mean (const long double w[], const size_t wstride, const long double data[], const size_t stride, const size_t n, const double mean);
double gsl_stats_long_double_wabsdev (const long double w[], const size_t wstride, const long double data[], const size_t stride, const size_t n);
double gsl_stats_long_double_wskew (const long double w[], const size_t wstride, const long double data[], const size_t stride, const size_t n);
double gsl_stats_long_double_wkurtosis (const long double w[], const size_t wstride, const long double data[], const size_t stride, const size_t n);

double gsl_stats_long_double_wcovariance (const long double w[], const size_t wstride, const long double data1[], const size_t stride1,const long double data2[], const size_t stride2, const size_t n);

double gsl_stats_long_double_wvariance_m (const long double w[], const size_t wstride, const long double data[], const size_t stride, const size_t n, const double wmean);
double gsl_stats_long_double_wsd_m (const long double w[], const size_t wstride, const long double data[], const size_t stride, const size_t n, const double wmean);
double gsl_stats_long_double_wabsdev_m (const long double w[], const size_t wstride, const long double data[], const size_t stride, const size_t n, const double wmean);
double gsl_stats_long_double_wskew_m_sd (const long double w[], const size_t wstride, const long double data[], const size_t stride, const size_t n, const double wmean, const double wsd);
double gsl_stats_long_double_wkurtosis_m_sd (const long double w[], const size_t wstride, const long double data[], const size_t stride, const size_t n, const double wmean, const double wsd);

double gsl_stats_long_double_wcovariance_m (const long double w[], const size_t wstride, const long double data1[], const size_t stride1,const long double data2[], const size_t stride2, const size_t n, const double mean1, const double mean2);

/* END OF FLOATING POINT TYPES */

double gsl_stats_long_double_pvariance (const long double data1[], const size_t stride1, const size_t n1, const long double data2[], const size_t stride2, const size_t n2);
double gsl_stats_long_double_ttest (const long double data1[], const size_t stride1, const size_t n1, const long double data2[], const size_t stride2, const size_t n2);

long double gsl_stats_long_double_max (const long double data[], const size_t stride, const size_t n);
long double gsl_stats_long_double_min (const long double data[], const size_t stride, const size_t n);
void gsl_stats_long_double_minmax (long double * min, long double * max, const long double data[], const size_t stride, const size_t n);

size_t gsl_stats_long_double_max_index (const long double data[], const size_t stride, const size_t n);
size_t gsl_stats_long_double_min_index (const long double data[], const size_t stride, const size_t n);
void gsl_stats_long_double_minmax_index (size_t * min_index, size_t * max_index, const long double data[], const size_t stride, const size_t n);

double gsl_stats_long_double_median_from_sorted_data (const long double sorted_data[], const size_t stride, const size_t n) ;
double gsl_stats_long_double_quantile_from_sorted_data (const long double sorted_data[], const size_t stride, const size_t n, const double f) ;

__END_DECLS

#endif /* __GSL_STATISTICS_LONG_DOUBLE_H__ */

[-- Attachment #4: gsl_statistics_float.h --]
[-- Type: text/plain, Size: 6637 bytes --]

/* statistics/gsl_statistics_float.h
 * 
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Jim Davies, Brian Gough
 * 
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#ifndef __GSL_STATISTICS_FLOAT_H__
#define __GSL_STATISTICS_FLOAT_H__

#include <stddef.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

double gsl_stats_float_mean (const float data[], const size_t stride, const size_t n);
double gsl_stats_float_variance (const float data[], const size_t stride, const size_t n);
double gsl_stats_float_sd (const float data[], const size_t stride, const size_t n);
double gsl_stats_float_variance_with_fixed_mean (const float data[], const size_t stride, const size_t n, const double mean);
double gsl_stats_float_sd_with_fixed_mean (const float data[], const size_t stride, const size_t n, const double mean);
double gsl_stats_float_absdev (const float data[], const size_t stride, const size_t n);
double gsl_stats_float_skew (const float data[], const size_t stride, const size_t n);
double gsl_stats_float_kurtosis (const float data[], const size_t stride, const size_t n);
double gsl_stats_float_lag1_autocorrelation (const float data[], const size_t stride, const size_t n);

double gsl_stats_float_covariance (const float data1[], const size_t stride1,const float data2[], const size_t stride2, const size_t n);

double gsl_stats_float_variance_m (const float data[], const size_t stride, const size_t n, const double mean);
double gsl_stats_float_sd_m (const float data[], const size_t stride, const size_t n, const double mean);
double gsl_stats_float_absdev_m (const float data[], const size_t stride, const size_t n, const double mean);
double gsl_stats_float_skew_m_sd (const float data[], const size_t stride, const size_t n, const double mean, const double sd);
double gsl_stats_float_kurtosis_m_sd (const float data[], const size_t stride, const size_t n, const double mean, const double sd);
double gsl_stats_float_lag1_autocorrelation_m (const float data[], const size_t stride, const size_t n, const double mean);

double gsl_stats_float_covariance_m (const float data1[], const size_t stride1,const float data2[], const size_t stride2, const size_t n, const double mean1, const double mean2);

/* DEFINED FOR FLOATING POINT TYPES ONLY */

double gsl_stats_float_wmean (const float w[], const size_t wstride, const float data[], const size_t stride, const size_t n);
double gsl_stats_float_wvariance (const float w[], const size_t wstride, const float data[], const size_t stride, const size_t n);
double gsl_stats_float_wsd (const float w[], const size_t wstride, const float data[], const size_t stride, const size_t n);
double gsl_stats_float_wvariance_with_fixed_mean (const float w[], const size_t wstride, const float data[], const size_t stride, const size_t n, const double mean);
double gsl_stats_float_wsd_with_fixed_mean (const float w[], const size_t wstride, const float data[], const size_t stride, const size_t n, const double mean);
double gsl_stats_float_wabsdev (const float w[], const size_t wstride, const float data[], const size_t stride, const size_t n);
double gsl_stats_float_wskew (const float w[], const size_t wstride, const float data[], const size_t stride, const size_t n);
double gsl_stats_float_wkurtosis (const float w[], const size_t wstride, const float data[], const size_t stride, const size_t n);

double gsl_stats_float_wcovariance (const float w[], const size_t wstride, const float data1[], const size_t stride1,const float data2[], const size_t stride2, const size_t n);

double gsl_stats_float_wvariance_m (const float w[], const size_t wstride, const float data[], const size_t stride, const size_t n, const double wmean);
double gsl_stats_float_wsd_m (const float w[], const size_t wstride, const float data[], const size_t stride, const size_t n, const double wmean);
double gsl_stats_float_wabsdev_m (const float w[], const size_t wstride, const float data[], const size_t stride, const size_t n, const double wmean);
double gsl_stats_float_wskew_m_sd (const float w[], const size_t wstride, const float data[], const size_t stride, const size_t n, const double wmean, const double wsd);
double gsl_stats_float_wkurtosis_m_sd (const float w[], const size_t wstride, const float data[], const size_t stride, const size_t n, const double wmean, const double wsd);

double gsl_stats_float_wcovariance_m (const float w[], const size_t wstride, const float data1[], const size_t stride1,const float data2[], const size_t stride2, const size_t n, const double mean1, const double mean2);
/* END OF FLOATING POINT TYPES */

double gsl_stats_float_pvariance (const float data1[], const size_t stride1, const size_t n1, const float data2[], const size_t stride2, const size_t n2);
double gsl_stats_float_ttest (const float data1[], const size_t stride1, const size_t n1, const float data2[], const size_t stride2, const size_t n2);

float gsl_stats_float_max (const float data[], const size_t stride, const size_t n);
float gsl_stats_float_min (const float data[], const size_t stride, const size_t n);
void gsl_stats_float_minmax (float * min, float * max, const float data[], const size_t stride, const size_t n);

size_t gsl_stats_float_max_index (const float data[], const size_t stride, const size_t n);
size_t gsl_stats_float_min_index (const float data[], const size_t stride, const size_t n);
void gsl_stats_float_minmax_index (size_t * min_index, size_t * max_index, const float data[], const size_t stride, const size_t n);

double gsl_stats_float_median_from_sorted_data (const float sorted_data[], const size_t stride, const size_t n) ;
double gsl_stats_float_quantile_from_sorted_data (const float sorted_data[], const size_t stride, const size_t n, const double f) ;

__END_DECLS

#endif /* __GSL_STATISTICS_FLOAT_H__ */

[-- Attachment #5: wcovar_source.c --]
[-- Type: text/plain, Size: 4524 bytes --]

/* statistics/wcovar_source.c
 * 
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Jim Davies, Brian Gough
 * 
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#if 0
static double
FUNCTION(compute,wcovariance) (const BASE w[], const size_t wstride, const BASE data[], const size_t stride, const size_t n, const double wmean)
{
  /* takes a dataset and finds the weighted covariance */

  long double wcovariance = 0 ;
  long double W = 0;

  size_t i;

  /* find the sum of the squares */
  for (i = 0; i < n; i++)
    {
      BASE wi = w[i * wstride];

      if (wi > 0) {
        const long double delta = (data[i * stride] - wmean);
        W += wi ;
        wcovariance += (delta * delta - wcovariance) * (wi / W);
      }
    }

  return wcovariance ;
}
#endif

static double
FUNCTION(compute,wcovariance) (const BASE w[], const size_t wstride, 
                               const BASE data1[], const size_t stride1,
                               const BASE data2[], const size_t stride2,
                               const size_t n, 
                               const double mean1, const double mean2);

static double
FUNCTION(compute,factor) (const BASE w[], const size_t wstride, const size_t n);

static double
FUNCTION(compute,wcovariance) (const BASE w[], const size_t wstride, 
                               const BASE data1[], const size_t stride1,
                               const BASE data2[], const size_t stride2,
                               const size_t n, 
                               const double wmean1, const double wmean2)
{
  /* computes weighted covariance */
  /* internally computes the factor ``N/(N-1)'' which multiplies the weighted sum-of-squared difference */

  /* Q: originally this was long double. why? */
  /*    I changed it to BASE to improve performance on processors that don't really have long double, (eg. mine) */

  BASE wcovariance = 0.0 ;
  /* W accumulates the sum of w */
  BASE W = 0.0;        
  /* W2 accumulates sum of w^2 */
  BASE W2 = 0.0;       

  /* don't forget that return value is only double precision*/
  double rval; 

  size_t i;

  for (i = 0; i < n; i++)
    {
      const BASE wi = w[i * wstride];
      const BASE delta1 = data1[i * stride1] - wmean1;
      const BASE delta2 = data2[i * stride2] - wmean2;

      if (wi > 0)
        {
          W += wi ;
          W2 += wi * wi;
          wcovariance += (delta1 * delta2 - wcovariance) * (wi / W);
        }
    }

  /* return mean covariance, scaled by weighted "N/N-1" coefficient */

  /* Q: Dividing by zero is this function's way of saying you had no positive weights.
     Perhaps a gsl_error (GSL_EDOMAIN... ) might be in order?  */
  rval = (W*W > 0.0) ? wcovariance * (W*W) / (W*W - W2) : 0.0 ;
  return rval;
}

double 
FUNCTION(gsl_stats,wcovariance_m) (const BASE w[], const size_t wstride, 
                                   const BASE data1[], const size_t stride1,
                                   const BASE data2[], const size_t stride2,
                                   const size_t n,
                                   const double wmean1, const double wmean2)
{
  const double wcov = FUNCTION(compute,wcovariance) (w, wstride, data1, stride1, data2, stride2, n, wmean1, wmean2);
  return wcov;
}

double 
FUNCTION(gsl_stats,wcovariance) (const BASE w[], const size_t wstride, 
                                   const BASE data1[], const size_t stride1,
                                   const BASE data2[], const size_t stride2,
                                   const size_t n)
{
  const double wmean1 = FUNCTION(gsl_stats,wmean) (w, wstride, data1, stride1, n);
  const double wmean2 = FUNCTION(gsl_stats,wmean) (w, wstride, data2, stride2, n);
  const double wcov = FUNCTION(compute,wcovariance) (w, wstride, data1, stride1, data2, stride2, n, wmean1, wmean2);
  return wcov;
}

[-- Attachment #6: Makefile.am --]
[-- Type: text/plain, Size: 1261 bytes --]

## Process this file with automake to produce Makefile.in

noinst_LTLIBRARIES = libgslstatistics.la

pkginclude_HEADERS = gsl_statistics.h gsl_statistics_char.h gsl_statistics_double.h gsl_statistics_float.h gsl_statistics_int.h gsl_statistics_long.h gsl_statistics_long_double.h gsl_statistics_short.h gsl_statistics_uchar.h gsl_statistics_uint.h gsl_statistics_ulong.h gsl_statistics_ushort.h

INCLUDES= -I$(top_builddir) -I$(top_srcdir)

libgslstatistics_la_SOURCES =  mean.c variance.c absdev.c skew.c kurtosis.c lag1.c p_variance.c minmax.c ttest.c median.c covariance.c quantiles.c wmean.c wcovar.c wvariance.c wabsdev.c wskew.c wkurtosis.c

noinst_HEADERS = mean_source.c variance_source.c covariance_source.c absdev_source.c skew_source.c kurtosis_source.c lag1_source.c p_variance_source.c minmax_source.c ttest_source.c median_source.c quantiles_source.c wmean_source.c wcovar_source.c wvariance_source.c wabsdev_source.c wskew_source.c wkurtosis_source.c test_float_source.c test_int_source.c

check_PROGRAMS = test
TESTS = $(check_PROGRAMS)

test_SOURCES = test.c test_nist.c
test_LDADD = libgslstatistics.la ../sort/libgslsort.la ../ieee-utils/libgslieeeutils.la ../err/libgslerr.la ../test/libgsltest.la ../sys/libgslsys.la ../utils/libutils.la



[-- Attachment #7: test_wcov.c --]
[-- Type: text/plain, Size: 348 bytes --]

#include <gsl/gsl_statistics.h>

int main()
{
  double x[] = {10.4, 5.6, 3.1, 6.4, 21.7,
    3,4, 5, 6, 7,
    1,2,3,4,5
  };

  printf( "%lf\n", gsl_stats_wcovariance(x+10, 1, x, 1, x, 1, 5));
  printf( "%lf\n", gsl_stats_wcovariance(x+10, 1, x, 1, x+5, 1, 5));
  printf( "%lf\n", gsl_stats_wcovariance(x+10, 1, x+5, 1, x+5, 1, 5));
  return 0;
}

Re: gsl_stats_wcovariance

[Brian Gough]

James Bergstra writes:
 > I'd appreciate it if someone who knows about unbiased estimators of weighted
 > sample statistics would have a glance at my algorithm.  I haven't a clue-- I
 > implemented a straightforward extension of the wvariance routines.

I think we need to keep the function factor() to allow use of a fixed
mean -- the least work would be to make covariance a copy of wvariance*.c
with an extra set of arguments added for data2.

-- 
Brian Gough

Re: gsl_stats_wcovariance

[Brian Gough]

James Bergstra writes:
 > Please find attached (if the mailing list supports attachments) the files that I
 > modified.  They are from the statistics subdirectory of today's cvs.  In
 > addition to rewriting wcovar_source.c, I added wcovariance declarations to
 > headers (double, float, long double) and added wcovar.c to Makefile.am.  Also
 > attached is a very simple test program.  This program's output agrees with the
 > cov.wt routine from R.

That's good, we do need tests before adding anything.

-- 
Brian Gough