We use cookies to ensure you have the best browsing experience on our website. Please read our cookie policy for more information about how we use cookies.
C++ implementation. This has been a great learning process for me as I'm new to C++ but if you just wish to practice statistics using C++ and not spend two evenings on this problem, I recommend using the Matrix class, which was 90% of my time spent on this problem.
I calculated the matrix inverse using LU decompozition (which is faster than Gauss elimination method), it's the 'invertLU' method.
#include<cmath>#include<cstdio>#include<vector>#include<iostream>#include<algorithm>usingnamespacestd;classMatrix{public:Matrix(size_trows,size_tcols);double&operator()(size_ti,size_tj);doubleoperator()(size_ti,size_tj)const;Matrixtranspose();boolprintMatrix();vector<Matrix>decompLU();intdeterminant(size_tignoreCol,size_tignoreRow);MatrixmatrixMul(Matrixm);MatrixinvertLU(MatrixL,MatrixU);MatrixinvertLU();size_tget_mRows();size_tget_mCols();private:size_tmRows;size_tmCols;vector<double>mData;};Matrix::Matrix(size_trows,size_tcols):mRows(rows),mCols(cols),mData(rows*cols){}size_tMatrix::get_mRows(){returnmRows;}size_tMatrix::get_mCols(){returnmCols;}double&Matrix::operator()(size_ti,size_tj){returnmData.at(i*mCols+j);}doubleMatrix::operator()(size_ti,size_tj)const{returnmData[i*mCols+j];}MatrixMatrix::transpose(){MatrixtM(mCols,mRows);for(size_ti=0;i<mRows;i++){for(size_tj=0;j<mCols;j++){tM(j,i)=(*this)(i,j);}}returntM;}intdeterminant(size_tignoreCol,size_tignoreRow){return0;}vector<Matrix>Matrix::decompLU(){vector<Matrix>matrixes;Matrixl(mRows,mCols);Matrixu(mRows,mCols);if(mRows!=mCols){cout<<"Matrix is not square. Can't do LU decompozition"<<endl;returnmatrixes;}size_ti=0,j=0,k=0;for(i=0;i<mRows;i++){for(j=0;j<mCols;j++){if(j<i)l(j,i)=0;else{l(j,i)=(*this)(j,i);for(k=0;k<i;k++){l(j,i)=l(j,i)-l(j,k)*u(k,i);}}}for(j=0;j<mRows;j++){if(j<i)u(i,j)=0;elseif(j==i)u(i,j)=1;else{u(i,j)=(*this)(i,j)/l(i,i);for(k=0;k<i;k++){u(i,j)=u(i,j)-((l(i,k)*u(k,j))/l(i,i));}}}}matrixes.push_back(l);matrixes.push_back(u);returnmatrixes;}MatrixMatrix::matrixMul(Matrixm){doublesum=0;MatrixmulRes(mRows,m.get_mCols());for(size_ti=0;i<mRows;i++){for(size_tj=0;j<m.get_mCols();j++){for(size_tk=0;k<mCols;k++){sum+=(*this)(i,k)*m(k,j);}mulRes(i,j)=sum;sum=0;}}returnmulRes;}boolMatrix::printMatrix(){for(size_ti=0;i<mRows;i++){for(size_tj=0;j<mCols;j++){std::cout<<mData.at(i*mCols+j)<<" ";}std::cout<<std::endl;}returntrue;}structLU_matrix{MatrixL;MatrixU;};MatrixMatrix::invertLU(MatrixL,MatrixU){MatrixeigenMatrix(L.get_mRows(),L.get_mCols());for(size_ti=0;i<eigenMatrix.get_mRows();i++){for(size_tj=0;j<eigenMatrix.get_mCols();j++){if(i==j){eigenMatrix(i,j)=1;}else{eigenMatrix(i,j)=0;}}}// calculate [L][Z] = [I] -> get three vectorsMatrixZ(L.get_mRows(),L.get_mCols());doubles;for(size_tZ_col=0;Z_col<L.get_mCols();Z_col++){// Iterate for each element (row) in Zfor(size_trow=0;row<L.get_mRows();row++){s=0;// On each row, iterate through existing Z elementsfor(size_tj=0;j<row;j++){s+=Z(j,Z_col)*L(row,j);}Z(row,Z_col)=(eigenMatrix(row,Z_col)-s)/L(row,row);}}// calculate [U][X] = [Z] -> get three vectors X, X = inv(A)MatrixX(U.get_mRows(),U.get_mCols());for(size_tX_col=0;X_col<U.get_mCols();X_col++){// Iterate for each element (row) in Z, go from bottom upsize_trow=0;for(size_ttemp_row=0;temp_row<U.get_mRows();temp_row++){row=U.get_mRows()-temp_row-1;//for (size_t row = U.get_mRows() - 1; row >= 0; row--) {s=0;// On each row, iterate through existing Z elementsfor(size_tj=row;j<U.get_mRows();j++){s+=X(j,X_col)*U(row,j);}X(row,X_col)=(Z(row,X_col)-s)/U(row,row);}}returnX;}MatrixMatrix::invertLU(){vector<Matrix>LU=(*this).decompLU();autoL=LU.at(0);autoU=LU.at(1);MatrixeigenMatrix(L.get_mRows(),L.get_mCols());for(size_ti=0;i<eigenMatrix.get_mRows();i++){for(size_tj=0;j<eigenMatrix.get_mCols();j++){if(i==j){eigenMatrix(i,j)=1;}else{eigenMatrix(i,j)=0;}}}// calculate [L][Z] = [I] -> get three vectorsMatrixZ(L.get_mRows(),L.get_mCols());doubles;for(size_tZ_col=0;Z_col<L.get_mCols();Z_col++){// Iterate for each element (row) in Zfor(size_trow=0;row<L.get_mRows();row++){s=0;// On each row, iterate through existing Z elementsfor(size_tj=0;j<row;j++){s+=Z(j,Z_col)*L(row,j);}Z(row,Z_col)=(eigenMatrix(row,Z_col)-s)/L(row,row);}}// calculate [U][X] = [Z] -> get three vectors X, X = inv(A)MatrixX(U.get_mRows(),U.get_mCols());for(size_tX_col=0;X_col<U.get_mCols();X_col++){// Iterate for each element (row) in Z, go from bottom upsize_trow=0;for(size_ttemp_row=0;temp_row<U.get_mRows();temp_row++){row=U.get_mRows()-temp_row-1;//for (size_t row = U.get_mRows() - 1; row >= 0; row--) {s=0;// On each row, iterate through existing Z elementsfor(size_tj=row;j<U.get_mRows();j++){s+=X(j,X_col)*U(row,j);}X(row,X_col)=(Z(row,X_col)-s)/U(row,row);}}returnX;}intmain(){boolp_intercept=true;size_tm,n;cin>>n>>m;doubletmp;if(p_intercept)n++;MatrixX(m,n);MatrixY(m,1);//load matrixfor(size_ti=0;i<m;i++){for(size_tj=0;j<n;j++){if(p_intercept&&j==0){X(i,j)=1.0;}else{cin>>tmp;X(i,j)=tmp;}}cin>>tmp;Y(i,0)=tmp;}cin>>m;MatrixX_test(m,n);for(size_ti=0;i<m;i++){for(size_tj=0;j<n;j++){if(p_intercept&&j==0){X_test(i,j)=1.0;}else{cin>>tmp;X_test(i,j)=tmp;}}}autoregressionWeights=((X.transpose()).matrixMul(X).invertLU()).matrixMul(X.transpose()).matrixMul(Y);//cout << "Testing regression weights" << endl;//regressionWeights.printMatrix();autoregressionPredict=X_test.matrixMul(regressionWeights);//cout << "Testing fitting" << endl;regressionPredict.printMatrix();}
Day 9: Multiple Linear Regression
You are viewing a single comment's thread. Return to all comments →
C++ implementation. This has been a great learning process for me as I'm new to C++ but if you just wish to practice statistics using C++ and not spend two evenings on this problem, I recommend using the Matrix class, which was 90% of my time spent on this problem.
I calculated the matrix inverse using LU decompozition (which is faster than Gauss elimination method), it's the 'invertLU' method.