如何使用GSL实现Pearson相关系数?

4
我有两个浮点数向量x和y,我想计算它们之间的Pearson相关系数。由于我需要对大量数据进行操作(例如10百万个不同的向量x和20千个不同的向量y),因此我使用C++编程语言,更具体地说是GSL中的gsl_stats_correlation函数。
以下是我的C++代码:
#include <iostream>
#include <vector>
using namespace std;

#include <gsl/gsl_vector.h>
#include <gsl/gsl_statistics.h>

int main (int argc, char ** argv)
{
  vector<double> x, y;
  size_t n = 5;
  x.push_back(1.0); y.push_back(1.0);
  x.push_back(3.1); y.push_back(3.2);
  x.push_back(2.0); y.push_back(1.9);
  x.push_back(5.0); y.push_back(4.9);
  x.push_back(2.0); y.push_back(2.1);
  for(size_t i=0; i<n; ++i)
      printf ("x[%ld]=%.1f y[%ld]=%.1f\n", i, x[i], i, y[i]);

  gsl_vector_const_view gsl_x = gsl_vector_const_view_array( &x[0], x.size() );
  gsl_vector_const_view gsl_y = gsl_vector_const_view_array( &y[0], y.size() );

  double pearson = gsl_stats_correlation( (double*) gsl_x.vector.data, sizeof(double),
                                          (double*) gsl_y.vector.data, sizeof(double),
                                          n );
  printf ("Pearson correlation = %f\n", pearson);

  return 0;
}

它成功编译(gcc -Wall -g pearson.cpp -lstdc++ -lgsl -lgslcblas -o pearson),但当我运行它时,输出如下:

x[0]=1.0 y[0]=1.0
x[1]=3.1 y[1]=3.2
x[2]=2.0 y[2]=1.9
x[3]=5.0 y[3]=4.9
x[4]=2.0 y[4]=2.1
Pearson correlation = 1.000000

以下是R代码,显然结果不应该完全为1:

x <- c(1.0,3.1,2.0,5.0,2.0); y <-c(1.0,3.2,1.9,4.9,2.1)
cor(x, y, method="pearson")  # 0.99798

我错过了什么?
1个回答

3

更改这行代码:

  double pearson = gsl_stats_correlation( (double*) gsl_x.vector.data, sizeof(double),
                                          (double*) gsl_y.vector.data, sizeof(double),
                                          n );

to:

  double pearson = gsl_stats_correlation( (double*) gsl_x.vector.data, 1,
                                          (double*) gsl_y.vector.data, 1,
                                          n );

如果您想避免重复使用“魔数”1:

  const size_t stride = 1;
  double pearson = gsl_stats_correlation( (double*) gsl_x.vector.data, stride,
                                          (double*) gsl_y.vector.data, stride,
                                          n );

gsl_stats_correlation假定使用double,第二个和第四个参数是“步长”的数量,因此通过给它sizeof(double),它会跳过sizeof(double)*sizeof(double)字节。


我根本不需要(double*)转换,但如果你需要它(也许你使用的是与我不同的gsl版本),你应该考虑使用static_cast<double*> - Ben Hocking
1
谢谢,我想不到。但现在我看不出有任何情况需要给一个不同于神奇的“1”的步幅。那么为什么需要两个步幅参数呢? - tflutre
2
我猜有些人可能有数据,其中每三个元素(例如)是相关信息。 - Ben Hocking

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接