算法#03--具体解释最小二乘法原理和代码

最小二乘法原理

最小二乘法的目标:求误差的最小平方和,相应有两种:线性和非线性。线性最小二乘的解是closed-form(例如以下文),而非线性最小二乘没有closed-form,通经常使用迭代法求解(如高斯牛顿迭代法,本文不作介绍)。

【首先得到线性方程组】

1.概念

最小二乘法(又称最小平方法)是一种数学优化技术。它通过最小化误差的平方和寻找数据的最佳函数匹配。

利用最小二乘法能够简便地求得未知的数据,并使得这些求得的数据与实际数据之间误差的平方和为最小。

最小二乘法还可用于曲线拟合。

2.原理

函数原型:

已知:

(x0,y0)。(x1。y1)…(xi,yi)…(xn,yn)个点,n>=k。

偏差平方和:

偏差平方和最小值能够通过使偏导数等于零得到:

简化左边等式有:

写成矩阵形式:公式①

将这个范德蒙得矩阵化简后可得到:公式②

也就是说X*A=Y,那么A = (X’*X)-1*X’*Y,便得到了系数矩阵A,同一时候,我们也就得到了拟合曲线。

高斯消元法

【然后解线性方程组,即公式①】

1.概念

数学上,高斯消元法(或译:高斯消去法)(英语:Gaussian Elimination),是线性代数中的一个算法,可用来为线性方程组求解,求出矩阵的秩,以及求出可逆方阵的逆矩阵。当用于一个矩阵时。高斯消元法会产生出一个“行梯阵式”。

2.原理

3.伪代码

这个算法和上面谈到的有点不同。它由绝对值最大的部分開始做起。这样能够改善算法的稳定性。

本算法由左至右地计算。每作出下面三个步骤,才跳到下一列和下一行:

  • 定出i列的绝对值最大的一个非0的数,将第i行的值与该行交换,使得该行拥有该列的最大值;
  • 将i列的数字除以该数,使得i列i行的数成为1。
  • 第(i+1)行下面(包含第(j+1)行)全部元素都转化为0。

全部步骤完毕后,这个矩阵会变成一个行梯矩阵,再用代入法就能够求解该方程组。

 i = 1
 j = 1
 while (i ≤ m and j ≤ n) do
   Find pivot in column j, starting in row i    // 从第i行開始。找出第j列中的最大值(i、j值应保持不变)
   maxi = i
   for k = i+1 to m do
     if abs(A[k,j]) > abs(A[maxi,j]) then
       maxi = k   // 使用交换法找出最大值(绝对值最大)
     end if
   end for
   if A[maxi,j] ≠ 0 then  // 判定找到的绝对值最大值是否为零:若不为零就进行下面操作;若为零则说明该列第(i+1)行下面(包含第(i+1)行)均为零,不须要再处理,直接跳转至第(j+1)列第(i+1)行
     swap rows i and maxi, but do not change the value of i   // 将第i行与找到的最大值所在行做交换,保持i值不变(i值记录了本次操作的起始行)
     Now A[i,j] will contain the old value of A[maxi,j].
     divide each entry in row i by A[i,j]    // 将交换后的第i行归一化(第i行全部元素分别除以A[i,j])
     Now A[i,j] will have the value 1.
     for u = i+1 to m do    // 第j列中,第(i+1)行下面(包含第(i+1)行)全部元素都减去A[i,j],直到第j列的i+1行以後元素均為零
       subtract A[u,j] * row i from row u
       Now A[u,j] will be 0, since A[u,j] - A[i,j] * A[u,j] = A[u,j] - 1 * A[u,j] = 0.
     end for
     i = i + 1
   end if
   j = j + 1  // 第j列中。第(i+1)行下面(包含第(i+1)行)全部元素均为零。

移至第(j+1)列,从第(i+1)行開始反复上述步骤。

end while

代码

public class CurveFitting {
     ///<summary>
    ///最小二乘法拟合二元多次曲线
    ///比如y=ax+b
    ///当中MultiLine将返回a。b两个參数。

///a相应MultiLine[1]
    ///b相应MultiLine[0]
    ///</summary>
    ///<param name="arrX">已知点的x坐标集合</param>
    ///<param name="arrY">已知点的y坐标集合</param>
    ///<param name="length">已知点的个数</param>
    ///<param name="dimension">方程的最高次数</param>
    public static double[] MultiLine(double[] arrX, double[] arrY, int length, int dimension) {
        int n = dimension + 1;                 //dimension次方程须要求 dimension+1个 系数
        double[][] Guass = new double[n][n + 1];
        for (int i = 0; i < n; i++){ //求矩阵公式①
            int j;
            for (j = 0; j < n; j++){
                Guass[i][j] = SumArr(arrX, j + i, length);//公式①等号左边第一个矩阵,即Ax=b中的A
            }
            Guass[i][j] = SumArr(arrX, i, arrY, 1, length);//公式①等号右边的矩阵,即Ax=b中的b
        }        

        return ComputGauss(Guass, n);//高斯消元法
    }

    //求数组的元素的n次方的和,即矩阵A中的元素
    private static double SumArr(double[] arr, int n, int length) {
        double s = 0;
        for (int i = 0; i < length; i++){
            if (arr[i] != 0 || n != 0){
                s = s + Math.pow(arr[i], n);
            }
            else{
                s = s + 1;
            }
        }
        return s;
    }

    //求数组的元素的n次方的和,即矩阵b中的元素
    private static double SumArr(double[] arr1, int n1, double[] arr2, int n2, int length) {
        double s = 0;
        for (int i = 0; i < length; i++)
        {
            if ((arr1[i] != 0 || n1 != 0) && (arr2[i] != 0 || n2 != 0))
                s = s + Math.pow(arr1[i], n1) * Math.pow(arr2[i], n2);
            else
                s = s + 1;
        }
        return s;
    }

    //高斯消元法解线性方程组
    private static double[] ComputGauss(double[][] Guass, int n) {
        int i, j;
        int k, m;
        double temp;
        double max;
        double s;
        double[] x = new double[n];

        for (i = 0; i < n; i++) {
            x[i] = 0.0;//初始化
        }

        for (j = 0; j < n; j++) {
            max = 0;
            k = j;

            // 从第i行開始,找出第j列中的最大值(i、j值应保持不变)
            for (i = j; i < n; i++) {
                if (Math.abs(Guass[i][j]) > max){
                    max = Guass[i][j];// 使用交换法找出最大值(绝对值最大)
                    k = i;
                }
            }

            if (k != j) {
                //将第j行与找到的最大值所在行做交换。保持i值不变(j值记录了本次操作的起始行)
                for (m = j; m < n + 1; m++) {
                    temp = Guass[j][m];
                    Guass[j][m] = Guass[k][m];
                    Guass[k][m] = temp;
                }
            }

            if (max == 0) {
                // "此线性方程为神秘线性方程"
                return x;
            }

            // 第m列中,第(j+1)行下面(包含第(j+1)行)全部元素都减去Guass[j][m] * s / (Guass[j][j])
            //直到第m列的i+1行以後元素均为零
            for (i = j + 1; i < n; i++) {
                s = Guass[i][j];
                for (m = j; m < n + 1; m++) {
                    Guass[i][m] = Guass[i][m] - Guass[j][m] * s / (Guass[j][j]);
                }
            }
        }//结束for (j=0;j<n;j++)

        //回代过程(见公式4.1.5)
        for (i = n - 1; i >= 0; i--) {
            s = 0;
            for (j = i + 1; j < n; j++) {
                s = s + Guass[i][j] * x[j];
            }
            x[i] = (Guass[i][n] - s) / Guass[i][i];
        }

        return x;
    }//返回值是函数的系数

    public static void main(String[] args)  {
        double[] x = {0, 1, 2, 3, 4, 5, 6, 7};
        double[] y = {0, 1, 4, 9, 16, 25, 36, 49};
        double[] a = MultiLine(x, y, 8, 2);

        for(int i =0; i <a.length;i++){
            System.out.println(a[i]);
        }
    }
}

输出:

0.708333333333342

-0.37500000000000583

1.0416666666666674

取整就得到y=x^2。

时间: 2024-10-12 16:51:39

算法#03--具体解释最小二乘法原理和代码的相关文章

算法学习#03--详解最小二乘法原理和代码

最小二乘法原理 最小二乘法的目标:求误差的最小平方和,对应有两种:线性和非线性.线性最小二乘的解是closed-form(如下文),而非线性最小二乘没有closed-form,通常用迭代法求解(如高斯牛顿迭代法,本文不作介绍). [首先得到线性方程组] 1.概念 最小二乘法(又称最小平方法)是一种数学优化技术.它通过最小化误差的平方和寻找数据的最佳函数匹配. 利用最小二乘法可以简便地求得未知的数据,并使得这些求得的数据与实际数据之间误差的平方和为最小. 最小二乘法还可用于曲线拟合. 2.原理 函

十大排序算法(原理及代码实现细节)

本文参考一些书籍啊哈算法,数据结构与算法(清华大学),已经一些网上的博客 然后动图也是从网上偷来的(^_^),代码实现我尽量用大家容易懂的方式实现 数组居多,然后,桶排序(是别人代码,不过写的不完全正确后面会更新),都是学习嘛 有误的地方,还望各位指正,希望对你有帮助(其实很灵活的,在运用上),也不要这样就满足了 多多地运用,会使理解更深的. 按上面的顺序来吧 原理在代码里直接上动图吧 冒泡排序动图演示 冒泡排序代码实现 1 #include<iostream> 2 #include<c

DeepLearning tutorial(3)MLP多层感知机原理简介+代码详解

DeepLearning tutorial(3)MLP多层感知机原理简介+代码详解 @author:wepon @blog:http://blog.csdn.net/u012162613/article/details/43221829 本文介绍多层感知机算法,特别是详细解读其代码实现,基于python theano,代码来自:Multilayer Perceptron,如果你想详细了解多层感知机算法,可以参考:UFLDL教程,或者参考本文第一部分的算法简介. 经详细注释的代码:放在我的gith

DeepLearning tutorial(4)CNN卷积神经网络原理简介+代码详解

DeepLearning tutorial(4)CNN卷积神经网络原理简介+代码详解 @author:wepon @blog:http://blog.csdn.net/u012162613/article/details/43225445 本文介绍多层感知机算法,特别是详细解读其代码实现,基于python theano,代码来自:Convolutional Neural Networks (LeNet).经详细注释的代码和原始代码:放在我的github地址上,可下载. 一.CNN卷积神经网络原理

DeepLearning tutorial(1)Softmax回归原理简介+代码详解

DeepLearning tutorial(1)Softmax回归原理简介+代码详解 @author:wepon @blog:http://blog.csdn.net/u012162613/article/details/43157801 本文介绍Softmax回归算法,特别是详细解读其代码实现,基于python theano,代码来自:Classifying MNIST digits using Logistic Regression,参考UFLDL. 一.Softmax回归简介 关于算法的详

十大经典排序算法最强总结(含Java代码实现)

最近几天在研究排序算法,看了很多博客,发现网上有的文章中对排序算法解释的并不是很透彻,而且有很多代码都是错误的,例如有的文章中在"桶排序"算法中对每个桶进行排序直接使用了Collection.sort()函数,这样虽然能达到效果,但对于算法研究来讲是不可以的.所以我根据这几天看的文章,整理了一个较为完整的排序算法总结,本文中的所有算法均有JAVA实现,经本人调试无误后才发出,如有错误,请各位前辈指出. 0.排序算法说明 0.1 排序的定义 对一序列对象根据某个关键字进行排序. 0.2

flume原理及代码实现

转载标明出处:http://www.cnblogs.com/adealjason/p/6240122.html 最近想玩一下流计算,先看了flume的实现原理及源码 源码可以去apache 官网下载 下面整理下flume的原理及代码实现: flume是一个实时数据收集工具,hadoop的生态圈之一,主要用来在分布式环境下各服务器节点做数据收集,然后汇总到统一的数据存储平台,flume支持多种部署架构模式,单点agent部署,分层架构模式部署,如通过一个负载均衡agent将收集的数据分发到各个子a

潜在语义分析Latent semantic analysis note(LSA)原理及代码实现

文章参考:http://blog.sina.com.cn/s/blog_62a9902f0101cjl3.html Latent Semantic Analysis (LSA)也被叫做Latent Semantic Indexing(LSI),从字面上的意思理解就是通过分析文档去发现这些文档中潜在的意思和概念.假设每个词仅表示一个概念,并且每个概念仅仅被一个词所描述,LSA将非常简单(从词到概念存在一个简单的映射关系) 不幸的是,这个问题并没有如此简单,因为存在不同的词表示同一个意思(同义词),

Base64加密解密原理以及代码实现(VC++)

Base64加密解密原理以及代码实现 转自:http://blog.csdn.net/jacky_dai/article/details/4698461 1. Base64使用A--Z,a--z,0--9,+,/ 这64个字符.    2. 编码原理:将3个字节转换成4个字节( (3 X 8) = 24 = (4 X 6) )先读入3个字节,每读一个字节,左移8位,再右移四次,每次6位,这样就有4个字节了.    3. 解码原理:将4个字节转换成3个字节.先读入4个6位(用或运算),每次左移6位