多元线性回归最小二乘法及其应用
Cholesky分解求系数参考:
[1]冯天祥. 多元线性回归最小二乘法及其经济分析[J]. 经济师,2003,11:129.
还可以采用最小二乘法来估计参数:
算法设计也可以参考两种系数最终公式设计。
下面的Java代码由网友设计,采用第一种方法计算参数。
1 /** 2 * 多元线性回归分析 3 * 4 * @param x[m][n] 5 * 每一列存放m个自变量的观察值 6 * @param y[n] 7 * 存放随即变量y的n个观察值 8 * @param m 9 * 自变量的个数 10 * @param n 11 * 观察数据的组数 12 * @param a 13 * 返回回归系数a0,...,am 14 * @param dt[4] 15 * dt[0]偏差平方和q,dt[1] 平均标准偏差s dt[2]返回复相关系数r dt[3]返回回归平方和u 16 * @param v[m] 17 * 返回m个自变量的偏相关系数 18 */ 19 public static void sqt2(double[][] x, double[] y, int m, int n, double[] a, 20 double[] dt, double[] v) { 21 int i, j, k, mm; 22 double q, e, u, p, yy, s, r, pp; 23 double[] b = new double[(m + 1) * (m + 1)]; 24 mm = m + 1; 25 b[mm * mm - 1] = n; 26 for (j = 0; j <= m - 1; j++) { 27 p = 0.0; 28 for (i = 0; i <= n - 1; i++) 29 p = p + x[j][i]; 30 b[m * mm + j] = p; 31 b[j * mm + m] = p; 32 } 33 for (i = 0; i <= m - 1; i++) 34 for (j = i; j <= m - 1; j++) { 35 p = 0.0; 36 for (k = 0; k <= n - 1; k++) 37 p = p + x[i][k] * x[j][k]; 38 b[j * mm + i] = p; 39 b[i * mm + j] = p; 40 } 41 a[m] = 0.0; 42 for (i = 0; i <= n - 1; i++) 43 a[m] = a[m] + y[i]; 44 for (i = 0; i <= m - 1; i++) { 45 a[i] = 0.0; 46 for (j = 0; j <= n - 1; j++) 47 a[i] = a[i] + x[i][j] * y[j]; 48 } 49 chlk(b, mm, 1, a); 50 yy = 0.0; 51 for (i = 0; i <= n - 1; i++) 52 yy = yy + y[i] / n; 53 q = 0.0; 54 e = 0.0; 55 u = 0.0; 56 for (i = 0; i <= n - 1; i++) { 57 p = a[m]; 58 for (j = 0; j <= m - 1; j++) 59 p = p + a[j] * x[j][i]; 60 q = q + (y[i] - p) * (y[i] - p); 61 e = e + (y[i] - yy) * (y[i] - yy); 62 u = u + (yy - p) * (yy - p); 63 } 64 s = Math.sqrt(q / n); 65 r = Math.sqrt(1.0 - q / e); 66 for (j = 0; j <= m - 1; j++) { 67 p = 0.0; 68 for (i = 0; i <= n - 1; i++) { 69 pp = a[m]; 70 for (k = 0; k <= m - 1; k++) 71 if (k != j) 72 pp = pp + a[k] * x[k][i]; 73 p = p + (y[i] - pp) * (y[i] - pp); 74 } 75 v[j] = Math.sqrt(1.0 - q / p); 76 } 77 dt[0] = q; 78 dt[1] = s; 79 dt[2] = r; 80 dt[3] = u; 81 } 82 private static int chlk(double[] a, int n, int m, double[] d) { 83 int i, j, k, u, v; 84 if ((a[0] + 1.0 == 1.0) || (a[0] < 0.0)) { 85 System.out.println("fail\n"); 86 return (-2); 87 } 88 a[0] = Math.sqrt(a[0]); 89 for (j = 1; j <= n - 1; j++) 90 a[j] = a[j] / a[0]; 91 for (i = 1; i <= n - 1; i++) { 92 u = i * n + i; 93 for (j = 1; j <= i; j++) { 94 v = (j - 1) * n + i; 95 a[u] = a[u] - a[v] * a[v]; 96 } 97 if ((a[u] + 1.0 == 1.0) || (a[u] < 0.0)) { 98 System.out.println("fail\n"); 99 return (-2); 100 } 101 a[u] = Math.sqrt(a[u]); 102 if (i != (n - 1)) { 103 for (j = i + 1; j <= n - 1; j++) { 104 v = i * n + j; 105 for (k = 1; k <= i; k++) 106 a[v] = a[v] - a[(k - 1) * n + i] * a[(k - 1) * n + j]; 107 a[v] = a[v] / a[u]; 108 } 109 } 110 } 111 for (j = 0; j <= m - 1; j++) { 112 d[j] = d[j] / a[0]; 113 for (i = 1; i <= n - 1; i++) { 114 u = i * n + i; 115 v = i * m + j; 116 for (k = 1; k <= i; k++) 117 d[v] = d[v] - a[(k - 1) * n + i] * d[(k - 1) * m + j]; 118 d[v] = d[v] / a[u]; 119 } 120 } 121 for (j = 0; j <= m - 1; j++) { 122 u = (n - 1) * m + j; 123 d[u] = d[u] / a[n * n - 1]; 124 for (k = n - 1; k >= 1; k--) { 125 u = (k - 1) * m + j; 126 for (i = k; i <= n - 1; i++) { 127 v = (k - 1) * n + i; 128 d[u] = d[u] - a[v] * d[i * m + j]; 129 } 130 v = (k - 1) * n + k - 1; 131 d[u] = d[u] / a[v]; 132 } 133 } 134 return (2); 135 } 136 /** 137 * @param args 138 */ 139 public static void main(String[] args) { 140 // TODO Auto-generated method stub 141 /** 142 * 一元回归 143 */ 144 // int i; 145 // double[] dt=new double[6]; 146 // double[] a=new double[2]; 147 // double[] x={ 0.0,0.1,0.2,0.3,0.4,0.5, 148 // 0.6,0.7,0.8,0.9,1.0}; 149 // double[] y={ 2.75,2.84,2.965,3.01,3.20, 150 // 3.25,3.38,3.43,3.55,3.66,3.74}; 151 // SPT.SPT1(x,y,11,a,dt); 152 // System.out.println(""); 153 // System.out.println("a="+a[1]+" b="+a[0]); 154 // System.out.println("q="+dt[0]+" s="+dt[1]+" p="+dt[2]); 155 // System.out.println(" umax="+dt[3]+" umin="+dt[4]+" u="+dt[5]); 156 /** 157 * 多元回归 158 */ 159 int i; 160 double[] a = new double[4]; 161 double[] v = new double[3]; 162 double[] dt = new double[4]; 163 double[][] x = { { 1.1, 1.0, 1.2, 1.1, 0.9 }, 164 { 2.0, 2.0, 1.8, 1.9, 2.1 }, { 3.2, 3.2, 3.0, 2.9, 2.9 } }; 165 double[] y = { 10.1, 10.2, 10.0, 10.1, 10.0 }; 166 SPT.sqt2(x, y, 3, 5, a, dt, v); 167 for (i = 0; i <= 3; i++) 168 System.out.println("a(" + i + ")=" + a[i]); 169 System.out.println("q=" + dt[0] + " s=" + dt[1] + " r=" + dt[2]); 170 for (i = 0; i <= 2; i++) 171 System.out.println("v(" + i + ")=" + v[i]); 172 System.out.println("u=" + dt[3]); 173 }
下面的C++代码由网友提供,采用第二中方法计算系数。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 | #include<iostream> #include<fstream> #include<iomanip> using namespace std; void transpose( double **p1, double **p2, int m, int n); void multipl( double **p1, double **p2, double **p3, int m, int n, int p); void Inver( double **p1, double **p2, int n); double SD( double **p1, double **p2, double **p3, double **p4, int m, int n); double ST( double **p1, int m); void de_allocate( double **data, int m); int main() { int row,col; char filename[30]; double SDsum,STsum,F,R2; cout<< "Input original data file: \n" ; ifstream infile; //打开文件 cin>>filename; infile.open(filename); if (!infile) { cout<< "Opening the file failed!\n" ; exit (1); } infile>>row>>col; //读入文件中的行数和列数 double **matrix= new double *[row]; //为动态二维数组分配内存 double **X= new double *[row]; double **Y= new double *[row]; double **XT= new double *[col]; double **XTX= new double *[col]; double **XTXInv= new double *[col]; double **XTXInvXT= new double *[col]; double **B= new double *[col]; double **YE= new double *[row]; for ( int i=0;i<row;i++) { matrix[i]= new double [col]; X[i]= new double [col]; Y[i]= new double [1]; Y[i]= new double [1]; YE[i]= new double [1]; } for ( int i=0;i<col;i++) { XT[i]= new double [row]; XTX[i]= new double [2*col]; /////////////////////为什么必须分配2*col列空间而不是col?在矩阵求逆时,XTX变增广矩阵,列数变为原来2吧倍,跟求逆算法有关。 XTXInv[i]= new double [col]; XTXInvXT[i]= new double [row]; B[i]= new double [1]; } for ( int i=0;i<row;i++) for ( int j=0;j<col;j++) infile>>matrix[i][j]; infile.close(); for ( int i=0;i<row;i++) { //提取1X和Y数组列 X[i][0]=1; Y[i][0]=matrix[i][col-1]; for ( int j=0;j<col-1;j++) X[i][j+1]=matrix[i][j]; } transpose(X,XT,row,col); multipl(XT,X,XTX,col,row,col); Inver(XTX,XTXInv,col); multipl(XTXInv,XT,XTXInvXT,col,col,row); multipl(XTXInvXT,Y,B,col,row,1); SDsum=SD(Y,X,B,YE,row,col); STsum=ST(Y,row); F=((STsum-SDsum)/(col-1))/(SDsum/(row-col)); R2=1/(1+(row-col)/F/(col-1)); cout<< "输出B:\n" ; //屏幕输出结果B,SD,ST,F,R2 for ( int i=0;i<col;i++) cout<<setiosflags(ios::fixed)<<setprecision(4)<<B[i][0]<< ' ' ; cout<<endl; cout<< "SD=" <<SDsum<< ';' << "ST=" <<STsum<< ';' << "F=" <<F<< ';' << "R2=" <<R2<<endl; ofstream outfile; // 结果写入文件 cout<< "Output file'name:\n" ; cin>>filename; outfile.open(filename); if (!outfile) { cout<< "Opening the file failed!\n" ; exit (1); } outfile<< "输出B:\n" ; for ( int i=0;i<col;i++) outfile<<B[i][0]<< ' ' ; outfile<<endl; outfile<<setiosflags(ios::fixed)<<setprecision(4)<< "SD=" <<SDsum<< ';' << "ST=" <<STsum<< ';' << "F=" <<F<< ';' << "R2=" <<R2<<endl; outfile<< "Y and YE and Y-YE's value are:\n" ; for ( int i=0;i<row;i++) outfile<<Y[i][0]<< " " <<YE[i][0]<< " " <<Y[i][0]-YE[i][0]<<endl; outfile.close(); de_allocate(matrix,row); de_allocate(X,row); de_allocate(Y,row); de_allocate(XT,col); de_allocate(XTX,col); de_allocate(XTXInv,col); de_allocate(XTXInvXT,col); de_allocate(B,col); de_allocate(YE,row); system ( "pause" ); return (0); } void de_allocate( double **data, int m) { //释放内存单元 for ( int i=0;i<m;i++) delete []data[i]; delete []data; } double ST( double **p1, int m) { //求总离差平方和ST double sum1=0,sum2=0,Yave=0; for ( int i=0;i<m;i++) sum1+=p1[i][0]; Yave=sum1/m; for ( int i=0;i<m;i++) sum2+=(p1[i][0]-Yave)*(p1[i][0]-Yave); return sum2; } double SD( double **p1, double **p2, double **p3, double **p4, int m, int n) { //求偏差平方和SD double sum1=0,sum2=0; for ( int i=0;i<m;i++) { sum1=0; for ( int k=0;k<n;k++) sum1+=p2[i][k]*p3[k][0]; p4[i][0]=sum1; } for ( int i=0;i<m;i++) sum2+=(p1[i][0]-p4[i][0])*(p1[i][0]-p4[i][0]); return sum2; } void transpose( double **p1, double **p2, int m, int n) { //矩阵转置 for ( int i=0;i<n;i++) for ( int j=0;j<m;j++) p2[i][j]=p1[j][i]; } void multipl( double **p1, double **p2, double **p3, int m, int n, int p) { //矩阵相乘 double sum; for ( int i=0;i<m;i++) { for ( int j=0;j<p;j++) { sum=0; for ( int k=0;k<n;k++) sum+=p1[i][k]*p2[k][j]; p3[i][j]=sum; } } } void Inver( double **p1, double **p2, int n) { //求逆矩阵 //初始化矩阵在右侧加入单位阵 for ( int i=0;i<n;i++) { for ( int j=0;j<n;j++) { p1[i][j+n]=0; p1[i][i+n]=1; } } //对于对角元素为0的进行换行操作 for ( int i=0;i<n;i++) { while (p1[i][i]==0) { for ( int j=i+1;j<n;j++) { if (p1[j][i]!=0) { double temp=0; for ( int r=i;r<2*n;r++) {temp=p1[j][r];p1[j][r]=p1[i][r];p1[i][r]=temp;} } break ; } } //if (p1[i][i]==0) return 0; } //行变换为上三角矩阵 double k=0; for ( int i=0;i<n;i++) { for ( int j=i+1;j<n;j++) { k=(-1)*p1[j][i]/p1[i][i]; for ( int r=i;r<2*n;r++) p1[j][r]+=k*p1[i][r]; } } //行变换为下三角矩阵 //double k=0; for ( int i=n-1;i>=0;i--) { for ( int j=i-1;j>=0;j--) { k=(-1)*p1[j][i]/p1[i][i]; for ( int r=0;r<2*n;r++) p1[j][r]+=k*p1[i][r]; } } //化为单位阵 for ( int i=n-1;i>=0;i--) { k=p1[i][i]; for ( int j=0;j<2*n;j++) p1[i][j]/=k; } //拆分出逆矩阵 for ( int i=0;i<n;i++) { for ( int j=0;j<n;j++) p2[i][j]=p1[i][n+j]; } } |
二元线性回归最小二乘法拟合空间直线。网友提供
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 | #include "stdafx.h" using namespace std; using std::cout; using std::cin; using std::endl; #include<math.h> #include <windows.h> /* 这是一个控制台程序,任何一个空间直线方程都能以如下的方式表示 x=az+b y=cz+d z=z 即 (x-b)/a=(y-d)/c=z/1 */ #include <vector> using std::vector; /* the fitting vaule of this line are: a=2 b=3 c=3 d=0 double z[10]={0.5,0.7,1,1.2,1.5,1.8,2,2.5,2.8,3}; double y[10]={1.5,2.1,3,3.6,4.5,5.4,6,7.5,8.4,9}; double x[10]={4,4.4,5,5.4,6,6.6,7,8,8.6,9}; */ double x[]={1,1.5,2,2.5,3,3.5,4,4.5,5}; double y[]={-8.1,-7.2,-6.2,-5.5,-4.8,-3.8,-3,-2.2,-1.3}; double z[]={-12,-11.8,-10.7,-9.5,-8.2,-7,-6,-4.5,-3.5}; int _tmain( int argc, _TCHAR* argv[]) { vector< double >position_z; vector< double >position_y; vector< double >position_x; if (( sizeof (z)/ sizeof ( double ))!=( sizeof (x)/ sizeof ( double ))||( sizeof (z)/ sizeof ( double ))!=( sizeof (y)/ sizeof ( double ))||( sizeof (x)/ sizeof ( double ))!=( sizeof (y)/ sizeof ( double ))) { ::MessageBox(NULL, "请检查输入数组的长度是否相等" , "方程无解!" ,MB_OK); } //int ss=sizeof(z)/sizeof(double); for ( int i=0;i<( sizeof (z)/ sizeof ( double ));++i) { position_z.push_back(z[i]); position_y.push_back(y[i]); position_x.push_back(x[i]); } //double m = position_z.size(); //caculate o1,p1,c1,o2,p2,c2 double o1 = 0.0; double p1 = 0.0; double x1 = 0.0; double o2 = 0.0; double p2 = 0.0; double x2 = 0.0; double y1 = 0.0; double y2 = 0.0; for (vector< double >::size_type t=0;t< position_z.size();++t) { o1 += position_z[t]*position_z[t]; p1 += position_z[t]; x1 += position_x[t]*position_z[t]; o2 += position_z[t]; p2 = position_z.size(); x2 += position_x[t] ; y2 += position_y[t]; y1 += position_y[t]*position_z[t]; } //caculate a b double a = 0.0 ,b = 0.0, c = 0.0, d = 0.0 ; if ((o1*p2-o2*p1)!= 0) { a = (x1*p2 - x2*p1) /(o1*p2-o2*p1); if (a<1e-10) { a=0; } b = (x2*o1 - x1*o2) /(o1*p2-o2*p1); if (b<1e-10) { b=0; } c = (y1*p2 - y2*p1) /(o1*p2-o2*p1); if (c<1e-10) { c=0; } d = (y2*o1 - y1*o2) /(o1*p2-o2*p1); if (d<1e-10) { d=0; } } else { ::MessageBox(NULL, "OK" , "方程无解!" ,MB_OK); } cout<<a<< " " <<b<< " " <<c<< " " <<d<<endl; system ( "pause" ); return 0; } |
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】凌霞软件回馈社区,博客园 & 1Panel & Halo 联合会员上线
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步