/**
* 多元线性回归分析
*
* @param x[m][n]
* 每一列存放m个自变量的观察值
* @param y[n]
* 存放随即变量y的n个观察值
* @param m
* 自变量的个数
* @param n
* 观察数据的组数
* @param a
* 返回回归系数a0,...,am
* @param dt[4]
* dt[0]偏差平方和q,dt[1] 平均标准偏差s dt[2]返回复相关系数r dt[3]返回回归平方和u
* @param v[m]
* 返回m个自变量的偏相关系数
*/
public static void sqt2(double[][] x, double[] y, int m, int n, double[] a,
double[] dt, double[] v) {
int i, j, k, mm;
double q, e, u, p, yy, s, r, pp;
double[] b = new double[(m + 1) * (m + 1)];
mm = m + 1;
b[mm * mm - 1] = n;
for (j = 0; j <= m - 1; j++) {
p = 0.0;
for (i = 0; i <= n - 1; i++)
p = p + x[j];
b[m * mm + j] = p;
b[j * mm + m] = p;
}
for (i = 0; i <= m - 1; i++)
for (j = i; j <= m - 1; j++) {
p = 0.0;
for (k = 0; k <= n - 1; k++)
p = p + x[k] * x[j][k];
b[j * mm + i] = p;
b[i * mm + j] = p;
}
a[m] = 0.0;
for (i = 0; i <= n - 1; i++)
a[m] = a[m] + y;
for (i = 0; i <= m - 1; i++) {
a = 0.0;
for (j = 0; j <= n - 1; j++)
a = a + x[j] * y[j];
}
chlk(b, mm, 1, a);
yy = 0.0;
for (i = 0; i <= n - 1; i++)
yy = yy + y / n;
q = 0.0;
e = 0.0;
u = 0.0;
for (i = 0; i <= n - 1; i++) {
p = a[m];
for (j = 0; j <= m - 1; j++)
p = p + a[j] * x[j];
q = q + (y - p) * (y - p);
e = e + (y - yy) * (y - yy);
u = u + (yy - p) * (yy - p);
}
s = Math.sqrt(q / n);
r = Math.sqrt(1.0 - q / e);
for (j = 0; j <= m - 1; j++) {
p = 0.0;
for (i = 0; i <= n - 1; i++) {
pp = a[m];
for (k = 0; k <= m - 1; k++)
if (k != j)
pp = pp + a[k] * x[k];
p = p + (y - pp) * (y - pp);
}
v[j] = Math.sqrt(1.0 - q / p);
}
dt[0] = q;
dt[1] = s;
dt[2] = r;
dt[3] = u;
}
private static int chlk(double[] a, int n, int m, double[] d) {
int i, j, k, u, v;
if ((a[0] + 1.0 == 1.0) || (a[0] < 0.0)) {
System.out.println("fail\n");
return (-2);
}
a[0] = Math.sqrt(a[0]);
for (j = 1; j <= n - 1; j++)
a[j] = a[j] / a[0];
for (i = 1; i <= n - 1; i++) {
u = i * n + i;
for (j = 1; j <= i; j++) {
v = (j - 1) * n + i;
a = a - a[v] * a[v];
}
if ((a + 1.0 == 1.0) || (a < 0.0)) {
System.out.println("fail\n");
return (-2);
}
a = Math.sqrt(a);
if (i != (n - 1)) {
for (j = i + 1; j <= n - 1; j++) {
v = i * n + j;
for (k = 1; k <= i; k++)
a[v] = a[v] - a[(k - 1) * n + i] * a[(k - 1) * n + j];
a[v] = a[v] / a;
}
}
}
for (j = 0; j <= m - 1; j++) {
d[j] = d[j] / a[0];
for (i = 1; i <= n - 1; i++) {
u = i * n + i;
v = i * m + j;
for (k = 1; k <= i; k++)
d[v] = d[v] - a[(k - 1) * n + i] * d[(k - 1) * m + j];
d[v] = d[v] / a;
}
}
for (j = 0; j <= m - 1; j++) {
u = (n - 1) * m + j;
d = d / a[n * n - 1];
for (k = n - 1; k >= 1; k--) {
u = (k - 1) * m + j;
for (i = k; i <= n - 1; i++) {
v = (k - 1) * n + i;
d = d - a[v] * d[i * m + j];
}
v = (k - 1) * n + k - 1;
d = d / a[v];
}
}
return (2);
}