最小二乘法拟合

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;
}