如何解决GSL中的多元线性回归正确吗?
我正在尝试使用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。
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");
////////////////////
}
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中,并执行多元线性回归。结果是正确的。这证明输入数据已正确读取。
解决方法
现在为gsl_multifit_linear
编写C驱动程序对我来说有点麻烦。据我所知,您必须从C或C ++运行它,因为R的gsl
程序包似乎不包含此函数的包装器...
但是,以您建议的两种方式(在R中)进行线性代数运算,得出的答案与lm()
相同:
m1 <- lm(Job_Perf~Mech_Apt+Consc,data=dd)
y <- dd$Job_Perf
X <- model.matrix(~Mech_Apt+Consc,data=dd)
cbind(lm=coef(m1),solve1=c(solve(crossprod(X)) %*% t(X) %*% y),solve2=c(solve(crossprod(X),t(X)%*%y)))
## lm solve1 solve2
## (Intercept) -4.10358123 -4.10358123 -4.10358123
## Mech_Apt 0.08640901 0.08640901 0.08640901
## Consc 0.08760164 0.08760164 0.08760164
在手动执行回归时,我看到的最常见错误是忘记R在模型矩阵中自动包含一个拦截列,但是当您列出一个alpha
估算值时,我想您并没有这样做错误...?
我最好的猜测是您以某种方式错误地设置了X
:看起来应该像这样
(Intercept) Mech_Apt Consc
1 40 25
1 45 20
1 38 30
...
数据:
dd <- read.table(header=TRUE,text="
Job_Perf Mech_Apt Consc
1 40 25
2 45 20
1 38 30
3 50 30
2 48 28
3 55 30
3 53 34
4 55 36
4 58 32
3 40 34
5 55 38
3 48 28
3 45 30
2 55 36
4 60 34
5 60 38
5 60 42
5 65 38
4 50 34
3 58 38
")
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。