多元线性回归最小二乘法及其应用

 
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 }
View Code
复制代码

下面的C++代码由网友提供,采用第二中方法计算系数。

 

二元线性回归最小二乘法拟合空间直线。网友提供

  



posted @   FeMiner  阅读(13135)  评论(0编辑  收藏  举报
点击右上角即可分享
微信分享提示