数值计算——最小二乘拟合二元一次多项式

￡神魔★判官ぃ 2022-06-17 03:53 251阅读 0赞

# 数值计算——最小二乘拟合二元一次多项式 #

## 最小二乘拟合： ##

就是根据一系列给定的数据点，求一条曲线使得数据点到曲线的某些（水平、竖直、垂直）距离最短。

## 推导过程： ##

1. 设拟合多项式为：

![1335515011_8365.jpg][]

2. 各点到这条曲线的距离之和，即偏差平方和如下：

![1335515269_4761.jpg][]

3. 为了求得符合条件的a值，对等式右边求ai偏导数，因而我们得到了：

![1335515297_9174.jpg][]

![1335515288_1734.jpg][]

.......

![1335515306_6536.jpg][]

4. 将等式左边进行一下化简，然后应该可以得到下面的等式：

![1335515389_2763.jpg][]

![1335515401_2536.jpg][]

.......

![1335515424_7676.jpg][]

5. 把这些等式表示成矩阵的形式，就可以得到下面的矩阵：

![1335515569_2057.jpg][]

6. 将这个范德蒙得矩阵化简后可得到:

![1335515598_9211.jpg][]

## 实现代码： ##

（里面用到的一些方法是前几篇博客写的）

下面这个代码是拟合y=ax+b的，拟合条件是数据点与拟合曲线的竖直距离最短，如果想要拟合水平距离最短的话只需要将原式写成x=cy+d同样的方法即可。也可以使用SVD方法，参见《科学计算导论》第3章。

package com.kexin.lab5;
    import com.kexin.lab4.NomalEquation;
    
    /**
     * 数据点与拟合直线竖直距离最短
     * @author KeXin
     *
     */
    public class LeastSquarel1 {
    	/**
    	 * 生成矩阵A
    	 * @return
    	 */
    	public static double[][] MakeA(double[] x){
    		double[][] a = new double[12][2];
    		for(int i = 0;i<12;i++){
    			for(int j = 0;j<2;j++){
    				a[i][j] = Math.round(Math.pow(x[i], j));
    			}
    		}
    		return a;
    	}
    	public static void main(String[] args) {
    		double[] x1= {0.5,1.1,1.7,2.1,2.5,2.9,3.3,3.7,4.2,4.9,5.3,6.0};//数据点x
    		double[] b1 = {1.6,2.4,3.8,4.3,4.7,4.8,5.5,6.1,6.3,7.1,7.1,8.2};//数据点y
    		double[][] A = MakeA(x1);//系数矩阵
    		//正规方程组方法解方程Ax = b
    		double[][] a2 = NomalEquation.Transposition(A); // A的转置矩阵
    		// 输入方程矩阵b
    		double[][] a = NomalEquation.MultiEquation(a2, A); // 系数矩阵的转置和系数矩阵的乘积n*n
    		// 求解变化后的结果矩阵a2*b1
    		int n = a2.length;
    		int m = b1.length;
    		double[] b = new double[n];
    		for (int i = 0; i < n; i++) {
    			for (int j = 0; j < m; j++) {
    				b[i] += a2[i][j] * b1[j];
    			}
    		}
    		// 楚列斯基分解系数矩阵的转置和系数矩阵的乘积a
    		a = NomalEquation.Cholesky(a);
    		m = b.length;
    		n = a.length;
    		double x[] = new double[m];
    		// 前代计算
    		for (int j = 0; j < m; j++) {
    			if (a[j][j] != 0) {
    				x[j] = b[j] / a[j][j];
    				b[j] = x[j];
    			}
    			for (int i = j + 1; i < m; i++) {
    				b[i] = b[i] - a[i][j] * x[j];
    			}
    		}
    		// 将a进行转置
    		a = NomalEquation.Transposition(a);
    		// 回代计算
    		for (int j = m - 1; j >= 0; j--) {
    			if (a[j][j] != 0) {
    				x[j] = b[j] / a[j][j];
    			}
    			for (int i = 0; i <= j - 1; i++) {
    				b[i] = b[i] - a[i][j] * x[j];
    			}
    		}
    		// 输出结果
    //		NomalEquation.PrintArray("x", x);
    		System.out.println("a :"+x[1]);
    		System.out.println("b :"+x[0]);
    		double r=0;
    		for(int i = 0;i<x1.length;i++){
    			r+=Math.pow((x[0]+x[1]*x1[i]-b1[i]),2);
    		}
    		System.out.println("||r||:"+r);
    	}
    }

[1335515011_8365.jpg]: https://img-my.csdn.net/uploads/201204/27/1335515011_8365.jpg
[1335515269_4761.jpg]: https://img-my.csdn.net/uploads/201204/27/1335515269_4761.jpg
[1335515297_9174.jpg]: https://img-my.csdn.net/uploads/201204/27/1335515297_9174.jpg
[1335515288_1734.jpg]: https://img-my.csdn.net/uploads/201204/27/1335515288_1734.jpg
[1335515306_6536.jpg]: https://img-my.csdn.net/uploads/201204/27/1335515306_6536.jpg
[1335515389_2763.jpg]: https://img-my.csdn.net/uploads/201204/27/1335515389_2763.jpg
[1335515401_2536.jpg]: https://img-my.csdn.net/uploads/201204/27/1335515401_2536.jpg
[1335515424_7676.jpg]: https://img-my.csdn.net/uploads/201204/27/1335515424_7676.jpg
[1335515569_2057.jpg]: https://img-my.csdn.net/uploads/201204/27/1335515569_2057.jpg
[1335515598_9211.jpg]: https://img-my.csdn.net/uploads/201204/27/1335515598_9211.jpg