zoukankan      html  css  js  c++  java
  • 工控 程序使用最小二乘法

    http://www.cnblogs.com/gnuhpc/archive/2012/12/09/2810000.html

    【算法研究与实现】最小二乘法直线拟合

     

    作者:gnuhpc 
    出处:http://www.cnblogs.com/gnuhpc/

    1.原理 
    在现实中经常遇到这样的问题,一个函数并不是以某个数学表达式的形式给出,而是以一些自变量与因变量的对应表给出,老师讲课的时候举的个例子是犯罪人的身高和留下的脚印长,可以测出一些人的数据然后得到一张表,它反应的是一个函数,回归的意思就是将它还原成数学表达式,这个式子也称为经验表达式,之所以叫经验就是说它不完全是实际中的那样准确,是有一定偏差的,只是偏差很小罢了。

    最小二乘法 
        设经验 
    方程是y=F(x),方程中含有一些待定系数an,给出真实值{(xi,yi)|i=1,2,...n},将这些x,y值代入方程然后作 
    差,可以描述误差:yi-F(xi),为了考虑整体的误差,可以取平方和,之所以要平方是考虑到误差可正可负直接相加可以相互抵消,所以记误差为:

    e=∑(yi-F(xi))^2

        它是一个多元函数,有an共n个未知量,现在要求的是最小值。所以必然满足对各变量的偏导等于0,于是得到n个方程:

    de/da1=0 
    de/da2=0 
    ... 
    de/dan=0

    n个方程确定n个未知量为常量是理论上可以解出来的。用这种误差分析的方法进行回归方程的方法就是最小二乘法。

    线性回归 
    如果经验方程是线性的,形如y=ax+b,就是线性回归。按上面的分析,误差函数为:

    e=∑(yi-axi-b)^2

    各偏导为:

    de/da=2∑(yi-axi-b)xi=0 
    de/db=-2∑(yi-axi-b)=0

    于是得到关于a,b的线性方程组:

    (∑xi^2)a+(∑xi)b=∑yixi 
    (∑xi)a+nb=∑yi

    设A=∑xi^2,B=∑xi,C=∑yixi,D=∑yi,则方程化为:

    Aa+Bb=C 
    Ba+nb=D

    解出a,b得:

    a=(Cn-BD)/(An-BB) 
    b=(AD-CB)/(An-BB)
     
    这就是我们要进行的算法。 
    2.C++实现 
    /* 
    * ===================================================================================== 

    *       Filename:  nihe.cpp 

    *    Description:  A least square method for fitting a curve 

    *        Version:  1.0 
    *        Created:  03/21/2009 12:32:56 PM 
    *       Revision:  none 
    *       Compiler:  gcc 

    *         Author:  Futuredaemon (BUPT), gnuhpc@gmail.com 
    *        Company:  BUPT_UNITED 

    * ===================================================================================== 
    */ 
    #include  <stdlib.h> 
    #include  <iostream> 
    #include  <valarray> 
    using namespace std; 
    int main(int argc, char *argv[]) 

        int num = 0; 
        cout << " Input how many numbers you want to calculate:"; 
        cin >> num; 
        valarray<double> data_x(num); 
        valarray<double> data_y(num); 
        while( num ) 
        { 
            cout << "Input the "<< num <<" of x:"; 
            cin >> data_x[num-1]; 
            cout << "Input the "<< num <<" of y:"; 
            cin >> data_y[num-1]; 
            num--; 
        } 
        double A =0.0; 
        double B =0.0; 
        double C =0.0; 
        double D =0.0; 
        A = (data_x*data_x).sum(); 
        B = data_x.sum(); 
        C = (data_x*data_y).sum(); 
        D = data_y.sum(); 
        double k,b,tmp =0; 
        if(tmp=(A*data_x.size()-B*B)) 
        { 
            k = (C*data_x.size()-B*D)/tmp; 
            b = (A*D-C*B)/tmp; 
        } 
        else 
        { 
            k=1; 
            b=0; 
        } 
        cout <<"k="<<k<<endl; 
        cout <<"b="<<b<<endl; 
        return 0; 

    3.OpenCV结构实现 
    #include "cv.h" 
    #include <iostream> 
    using namespace std; 
    int main(int argc, char *argv[]) 

      int i=0; 
      int j=0; 
      int num; 
      double A,B,C,D; 
      double k,b,tmp=0; 
      cout <<"Input how many numbers you want to calculate:"; 
      cin >>num; 
      CvMat *mat1=cvCreateMat(1,num,CV_64FC1); 
      CvMat *mat2=cvCreateMat(1,num,CV_64FC1); 
      CvMat *mattmp=cvCreateMat(1,num,CV_64FC1); 
      for (j=0;j<mat1->cols;j++) 
        { 
          cout << "data X"<<j<<"="; 
          cin>>CV_MAT_ELEM(*mat1,double,0,j); 
          cout << "data Y"<<j<<"="; 
          cin>>CV_MAT_ELEM(*mat2,double,0,j); 
        } 
      for (j=0;j<mat1->cols;j++) 
        { 
          cout<<"X="<<CV_MAT_ELEM(*mat1,double,0,j) 
              <<",Y="<<CV_MAT_ELEM(*mat2,double,0,j)<<endl; 
        } 
      cvMul(mat1,mat1,mattmp,1); 
      A = cvSum(mattmp).val[0]; 
      B = cvSum(mat1).val[0]; 
      cvMul(mat1,mat2,mattmp,1); 
      C = cvSum(mattmp).val[0]; 
      D = cvSum(mat2).val[0]; 
      tmp = A*mat1->cols-B*B; 
      k = (C*mat1->cols-B*D)/tmp; 
      b = (A*D-C*B)/tmp; 
      cout << "k=" << k <<endl; 
      cout << "b=" << b <<endl; 
      cvReleaseMat(&mat1); 
      cvReleaseMat(&mat2); 
      return 0; 
    }

    作者:gnuhpc 
    出处:http://www.cnblogs.com/gnuhpc/

  • 相关阅读:
    Atitit 趋势管理之道 attilax著
    Atitit 循环处理的新特性 for...else...
    Atitit 2017年的技术趋势与未来的大技术趋势
    atitit 用什么样的维度看问题.docx 如何了解 看待xxx
    atitit prj mnrs 项目中的几种经理角色.docx
    Atitit IT办公场所以及度假村以及网点以及租房点建设之道 attilax总结
    Atitit 工具选型的因素与方法 attilax总结
    Atitit.团队文化建设影响组织的的一些原理 法则 定理 效应 p826.v4
    Atiitt 管理方面的误区总结 attilax总结
    Atitit 未来趋势把控的书籍 attilax总结 v3
  • 原文地址:https://www.cnblogs.com/zhangzhifeng/p/5980705.html
Copyright © 2011-2022 走看看