【发布时间】：2020-10-26 02:35:21
【问题描述】：
我尝试使用 GSL 执行多元线性回归。 beta 的估计量是 beta=(X'X)^(-1)X'y。使用了三种方法 -
- 直接，即使用所需操作计算此公式的右侧；
- 求解X'Xb=X'y，即求解系统Ax=b，其中A=X'X，b=X'y，x=b；
- 使用函数 gsl_multifit_linear。测试数据来自here。
所有三种方法都给出了相同的错误结果：
beta1=0.0365505，beta2=0.0435827，alpha=0.645627。
相同的数据文件在 R 中产生正确的输出：
beta1=0.086409，beta2=0.087602，alpha=-4.1。
另一个数据集的结果（来自here 的 WHO 预期寿命数据）在 GSL 和 R 之间也存在差异。
很可能我在某个地方有错误。但是，如果至少有人可以确认 GSL 多元线性回归正常工作，那就太好了。
编辑1。 语言是 C++。 option1的代码，即beta=(X'X)^(-1)X'y
void mols(vector <double> predictandVector,vector <vector <double> > predictorMatrix) {
//the function accepts vector of predictant values and matrix of predictands
//convert input vectors into gsl_vector
int n=predictandVector.size();
int m=predictorMatrix[0].size();
gsl_vector*y=gsl_vector_alloc(n);
gsl_matrix*x=gsl_matrix_alloc(n,m);
for(int i=0;i<n;i++) {
gsl_vector_set(y,i,predictandVector[i]);
for(int j=0;j<m;j++) {
gsl_matrix_set(x,i,j,predictorMatrix[i][j]);
}
}
//multiply X' by X
gsl_matrix*xxTranspose=gsl_matrix_alloc(m,m);
gsl_blas_dgemm(CblasTrans,CblasNoTrans,1.0,x,x,0.0,xxTranspose);
/////////////////////////////////
//perform LU decomposition of xxTranspose to find it's inverse
gsl_permutation*p=gsl_permutation_alloc(m);
int s;
gsl_linalg_LU_decomp(xxTranspose,p,&s);
//////////////////////////////
//compute inverse of xxTranspose matrix
gsl_matrix*xxTransposeInverse=gsl_matrix_alloc(m,m);
gsl_linalg_LU_invert(xxTranspose,p,xxTransposeInverse);
//////////////////////////////
//multiply inverse of xxTranspose matrix by xTranspose
gsl_matrix*xxTransposeInverseXtranspose=gsl_matrix_alloc(m,n);
gsl_blas_dgemm(CblasNoTrans,CblasTrans,1.0,xxTransposeInverse,x,0.0,xxTransposeInverseXtranspose);
//////////////////////////////////////
//multiply matrix (X'X)^(-1)X' and vector y; this must give values of beta
gsl_vector*b=gsl_vector_alloc(m);
gsl_blas_dgemv(CblasNoTrans,1.0,xxTransposeInverseXtranspose,y,0.0,b);
//write beta values to file
FILE*f;
f=fopen("molsResult.dat","w");
gsl_vector_fprintf(f,b,"%g");
////////////////////
}
option2的代码，即求解系统X'Xbeta=X'y
oid mols(vector <double> predictandVector,vector <vector <double> > predictorMatrix) {
//convert input vectors into gsl_vector
int n=predictandVector.size();
int m=predictorMatrix[0].size();
gsl_vector*y=gsl_vector_alloc(n);
gsl_matrix*x=gsl_matrix_alloc(n,m);
for(int i=0;i<n;i++) {
gsl_vector_set(y,i,predictandVector[i]);
for(int j=0;j<m;j++) {
gsl_matrix_set(x,i,j,predictorMatrix[i][j]);
}
}
/////////////////////
//compute X'X
gsl_matrix*xxTranspose=gsl_matrix_alloc(m,m);
gsl_blas_dgemm(CblasTrans,CblasNoTrans,1.0,x,x,0.0,xxTranspose);
/////////////////////////////////
//compute X'y
gsl_vector*XtransposeY=gsl_vector_alloc(m);
gsl_blas_dgemv(CblasTrans,1.0,x,y,0.0,XtransposeY);
//////////////////////////////////////
//solve the linear system X'Xb=X'y to find coefficients beta
gsl_vector*b=gsl_vector_alloc(m);
int k=n;
if(m<n) {
k=m;
}
gsl_vector*tau=gsl_vector_alloc(k);
gsl_linalg_QR_decomp(xxTranspose,tau);
gsl_linalg_QR_solve(xxTranspose,tau,XtransposeY,b);
//write beta values into file
FILE*f;
f=fopen("molsResult.dat","w");
gsl_vector_fprintf(f,b,"%g");
////////////////////////////////
}
option3的代码，即使用gsl_multifit_linear。
void mols(vector <double> predictandVector,vector <vector <double> > predictorMatrix) {
//convert input vectors into gsl_vector
int n=predictandVector.size();
int m=predictorMatrix[0].size();
gsl_vector*y=gsl_vector_alloc(n);
gsl_matrix*x=gsl_matrix_alloc(n,m);
for(int i=0;i<n;i++) {
gsl_vector_set(y,i,predictandVector[i]);
for(int j=0;j<m;j++) {
gsl_matrix_set(x,i,j,predictorMatrix[i][j]);
}
}
///////////////////////////
//Use gsl_multifit_linear
gsl_multifit_linear_workspace*w=gsl_multifit_linear_alloc(n,2);
double chisq;
gsl_matrix*cov=gsl_matrix_alloc(m,m);
gsl_vector*b=gsl_vector_alloc(m);
gsl_multifit_linear(x,y,b,cov,&chisq,w);
//////////////////////
//write beta values into file
FILE*f;
f=fopen("molsResult.dat","w");
gsl_vector_fprintf(f,b,"%g");
//////////////////////////
}
编辑2。 还要绝对确保正确读取输入数据，添加到 option3 代码中的以下代码：
FILE*xyWrite;
xyWrite=fopen("molsInputData.dat","w");
for(int i=0;i<n;i++) {
fputs(to_string(gsl_vector_get(y,i)).c_str(),xyWrite);
fputs(" ",xyWrite);
fputs(to_string(gsl_matrix_get(x,i,0)).c_str(),xyWrite);
fputs(" ",xyWrite);
fputs(to_string(gsl_matrix_get(x,i,1)).c_str(),xyWrite);
fputs("\n",xyWrite);
}
此代码将 gsl_vector(s) x 和 y（用作输入）写入文件。然后将该文件加载到 R 中并执行多元线性回归。结果是正确的。这证明输入数据被正确读入。
【问题讨论】：
-
请在问题中发布您尝试过的代码。 当您说 GSL 时，您是在运行 C 代码、C++ 代码还是 R 包
gsl？ -
也许this SO post 可以提供帮助。
标签： r linear-regression gsl