二元线性回归及java实现

建立多元线性回归方程,实际上是对多元线性模型(2-2-4)进行估计,寻求估计式(2-2-3)的过程。与一元线性回归(介绍见上篇)分析相同,其基本思想是根据最小二乘原理,求解 使全部观测值 与回归值 的残差平方和达到最小值。由于残差平方和

          (2-2-5)

    是 的非负二次式,所以它的最小值一定存在。

    根据极值原理,当Q取得极值时, 应满足

    由(2-2-5)式,即满足

                    (2-2-6)

    (2-2-6)式称为正规方程组。它可以化为以下形式

     (2-2-7)

    如果用A表示上述方程组的系数矩阵可以看出A是对称矩阵。则有

                   (2-2-8)

 

 式中X是多元线性回归模型中数据的结构矩阵, 是结构矩阵X的转置矩阵。

 (2-2-7)式右端常数项也可用矩阵D来表示

    即

           

    因此(2-2-7)式可写成

Ab=D           (2-2-10)

    或

            (2-2-11)

如果A满秩(即A的行列式  )那么A的逆矩阵A-1存在,则由(2-10)式和(2-11)式得 的最小二乘估计为

            (2-2-12)

  也就是多元线性回归方程的回归系数。

  为了计算方便往往并不先求  ,再求b,而是通过解线性方程组(2-2-7)来求b。(2-2-7)是一个有p+1个未知量的线性方程组,它的第一个方程可化为

            (2-2-13)

    式中

            (2-2-14)

    将(2-2-13)式代入(2-2-7)式中的其余各方程,得

            (2-2-15)

    其中

            (2-2-16)

    将方程组(2-2-15)式用矩阵表示,则有

Lb=F           (2-2-17)

    其中

  

    于是

b=L-1F           (2-2-18)

  因此求解多元线性回归方程的系数可由(2-2-16)式先求出L,然后将其代回(2-2-17)式中求解。求b时,可用克莱姆法则求解,也可通过高斯变换求解。如果把b直接代入(2-2-18)式,由于要先求出L的逆矩阵,因而相对复杂一些。

  例2-2-1 表2-2-1为某地区土壤内含植物可给态磷(y)与土壤内所含无机磷浓度(x1)、土壤内溶于K2CO3溶液并受溴化物水解的有机磷浓度(x2)以及土壤内溶于K2CO3溶液但不溶于溴化物的有机磷(x3)的观察数据。求y对x1, x2, x3的线性回归方程  。

表2-2-1 土壤含磷情况观察数据

    计算如下:

    由(2-2-16)式

  

  

 代入(2-2-15)式得

         (2-2-19)

    若用克莱姆法则解上述方程组,则其解为

                (2-2-20)

    其中

    计算得

b1=1.7848,b2=-0.0834,b3=0.1611

     回归方程为


java实现

/**
	 * 多元线性回归分析
	 * 
	 * @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][i];
			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[i][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[i];
		for (i = 0; i <= m - 1; i++) {
			a[i] = 0.0;
			for (j = 0; j <= n - 1; j++)
				a[i] = a[i] + x[i][j] * y[j];
		}
		chlk(b, mm, 1, a);
		yy = 0.0;
		for (i = 0; i <= n - 1; i++)
			yy = yy + y[i] / 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][i];
			q = q + (y[i] - p) * (y[i] - p);
			e = e + (y[i] - yy) * (y[i] - 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][i];
				p = p + (y[i] - pp) * (y[i] - 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[u] = a[u] - a[v] * a[v];
			}
			if ((a[u] + 1.0 == 1.0) || (a[u] < 0.0)) {
				System.out.println("fail\n");
				return (-2);
			}
			a[u] = Math.sqrt(a[u]);
			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[u];
				}
			}
		}
		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[u];
			}
		}
		for (j = 0; j <= m - 1; j++) {
			u = (n - 1) * m + j;
			d[u] = d[u] / 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[u] = d[u] - a[v] * d[i * m + j];
				}
				v = (k - 1) * n + k - 1;
				d[u] = d[u] / a[v];
			}
		}
		return (2);
	}

	/**
	 * @param args
	 */
	public static void main(String[] args) {
		// TODO Auto-generated method stub
		/**
		 * 一元回归
		 */
		// int i;
		// double[] dt=new double[6];
		// double[] a=new double[2];
		// double[] x={ 0.0,0.1,0.2,0.3,0.4,0.5,
		// 0.6,0.7,0.8,0.9,1.0};
		// double[] y={ 2.75,2.84,2.965,3.01,3.20,
		// 3.25,3.38,3.43,3.55,3.66,3.74};
		// SPT.SPT1(x,y,11,a,dt);
		// System.out.println("");
		// System.out.println("a="+a[1]+" b="+a[0]);
		// System.out.println("q="+dt[0]+" s="+dt[1]+" p="+dt[2]);
		// System.out.println(" umax="+dt[3]+" umin="+dt[4]+" u="+dt[5]);
		/**
		 * 多元回归
		 */
		int i;
		double[] a = new double[4];
		double[] v = new double[3];
		double[] dt = new double[4];

		double[][] x = { { 1.1, 1.0, 1.2, 1.1, 0.9 },
				{ 2.0, 2.0, 1.8, 1.9, 2.1 }, { 3.2, 3.2, 3.0, 2.9, 2.9 } };
		double[] y = { 10.1, 10.2, 10.0, 10.1, 10.0 };
		SPT.sqt2(x, y, 3, 5, a, dt, v);
		for (i = 0; i <= 3; i++)
			System.out.println("a(" + i + ")=" + a[i]);
		System.out.println("q=" + dt[0] + "  s=" + dt[1] + "  r=" + dt[2]);
		for (i = 0; i <= 2; i++)
			System.out.println("v(" + i + ")=" + v[i]);
		System.out.println("u=" + dt[3]);

	}


评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值