最小二乘法拟合

来源：互联网发布：在澳洲怎么用淘宝编辑：程序博客网时间：2024/04/28 00:16
main函数里面是拟合的数据。调用ercheng(i)输出i次的拟合曲线的系数。系数依次从高次到低次。
/***********************************************\ |Author: Messyidea |Created Time: 2014/4/29 20:26:02 |File Name: 最小二乘法.cpp |Description: \***********************************************/#include <iostream>#include <cstdio>#include <cmath>#include <cstdlib>#include <string>#include <cstring>#include <algorithm>#include <vector>#include <list>#include <map>#include <set>#include <deque>#include <queue>#include <stack>using namespace std;vector <pair<int,int> > v;//列主元消元 bool solve(double** data,int n,int m,double* ans){    double maxn;    int p;    for(int i=0;i<n;++i){        //找主元，并交换        maxn = data[i][i];        p = i;        for(int j=i+1;j<n;++j){            if(fabs(maxn) < fabs(data[j][i])){                p = j;                maxn = data[i][i];            }        }        if(i!=p)        {            for(int k=i;k<m;++k){                swap(data[i][k],data[p][k]);            }        }        //判断是否可解        if(data[i][i] == 0) return false;        //高斯消元        for(int j=i+1;j<n;++j){            double shang = data[j][i] / data[i][i];            for(int k=i;k<m;++k){                data[j][k] -= data[i][k]*shang;            }         }    }    //回代    int pos = 0;    for(int i=n-1;i>=0;--i){        double sum = data[i][m-1];        for(int j=0;j<pos;++j){            sum -= data[i][n-1-j] * ans[j];        }        sum/=data[i][i];        ans[pos++] = sum;    }    return true;}//次方 int ppow(int x,int n){    int sum = 1;    while(n -- ){        sum*=x;    }    return sum;}void ercheng(int num){    //初始化基础数据     int* sum_x = new int[num*2];    int* sum_xny = new int[num+1];    for(int i=0;i<num*2;++i){        sum_x[i] = 0;        for(int j=0;j<v.size();++j){            sum_x[i] += ppow(v[j].first,i+1);        }    }    //for(int i=0;i<num*2;++i) cout<<sum_x[i]<<" ";cout<<endl;    int** ma = new int*[num+1];    for(int i=0;i<num+1;++i) ma[i] = new int[num+1];    ma[0][0] = 9;    for(int i=1;i<num+1;++i) ma[0][i] = sum_x[i-1];    for(int i=1;i<num+1;++i){        for(int j=0;j<num+1;++j) ma[i][j] = sum_x[j+i-1];    }    for(int i=0;i<num+1;++i){        sum_xny[i] = 0;        for(int j=0;j<v.size();++j){            sum_xny[i] += ppow(v[j].first,i)*v[j].second;        }    }    //构造增广矩阵     double** data = new double*[num+1];    for(int i=0;i<num+1;++i) data[i] = new double[num+2];    for(int i=0;i<num+1;++i){        for(int j=0;j<num+1;++j){            data[i][j] = ma[i][j];        }    }    for(int i=0;i<num+1;++i) data[i][num+1] = sum_xny[i];    /*for(int i=0;i<num+1;++i){        for(int j=0;j<num+2;++j){            cout<<data[i][j]<<" ";        }        cout<<endl;    }*/    double* ans = new double[num+1];    //列主元消元法     solve(data,num+1,num+2,ans);    for(int i=0;i<num+1;++i){        cout<<ans[i]<<" ";    }    cout<<endl;}void init(){    //初始化数据     v.push_back(make_pair(1,10));    v.push_back(make_pair(3,5));    v.push_back(make_pair(4,4));    v.push_back(make_pair(5,2));    v.push_back(make_pair(6,1));    v.push_back(make_pair(7,1));    v.push_back(make_pair(8,2));    v.push_back(make_pair(9,3));    v.push_back(make_pair(10,4));}int main() {    //freopen("input.txt","r",stdin);     //初始化数据    init();    for(int i=1;i<8;++i){        cout<<"cishu "<<i<<endl;        ercheng(i);    }     return 0;}
0 0