GSL中的多元线性回归正确吗？

编程问答 2022-08-09

问题描述

我正在尝试使用GSL执行多元线性回归。 beta的估计值是beta =（X'X）^（-1）X'y。使用了三种方法-

  1. 直接，即使用所需的运算来计算该公式的右侧；
  2. 求解X'Xb = X'y，即求解系统Ax = b，其中A = X'X，b = X'y和x = b；
  3. 使用函数gsl_multifit_linear。测试数据来自here。

这三种方法都给出了错误的结果：

beta1 = 0.0365505，beta2 = 0.0435827，alpha = 0.645627。

同一数据文件在R中产生正确的输出：

beta1 = 0.086409，beta2 = 0.087602，alpha = -4.1。

GSL和R之间的另一个数据集（来自here的WHO预期寿命数据）的结果也有所不同。

我很可能在某个地方犯了一个错误。但是，如果至少有人可以确认GSL多元线性回归能够正常工作，那将是很好的。

编辑1。 语言是C ++。 选项1的代码，即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,j,predictorMatrix[i][j]);
                }
        }
        //multiply X' by X
        gsl_matrix*xxTranspose=gsl_matrix_alloc(m,m);
        gsl_blas_dgemm(CblasTrans,CblasNoTrans,1.0,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,xxTransposeInverse);
        //////////////////////////////
        //multiply inverse of xxTranspose matrix by xTranspose
        gsl_matrix*xxTransposeInverseXtranspose=gsl_matrix_alloc(m,n);
         gsl_blas_dgemm(CblasNoTrans,CblasTrans,xxTransposeInverse,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,xxTransposeInverseXtranspose,y,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,predictorMatrix[i][j]);
                }
        }
        /////////////////////
        //compute X'X
        gsl_matrix*xxTranspose=gsl_matrix_alloc(m,xxTranspose);
        /////////////////////////////////
        //compute X'y
        gsl_vector*XtransposeY=gsl_vector_alloc(m);
        gsl_blas_dgemv(CblasTrans,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","%g");
        ////////////////////////////////
}

option3的代码，即使用gsl_multifit_linear。

void mols(vector <double> predictandVector,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,cov,&chisq,w);
        //////////////////////
        //write beta values into file
        FILE*f;
        f=fopen("molsResult.dat","%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,0)).c_str(),1)).c_str(),xyWrite);
                fputs("\n",xyWrite);
        }

此代码将gsl_vector（s）x和y（用作输入）写入文件。然后将此文件加载到R中，并执行多元线性回归。结果是正确的。这证明输入数据已正确读取。

解决方法

暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！

如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。

小编邮箱:dio#foxmail.com (将#修改为@）

gsllinear-regressionrr