grm.cpp

多项式曲线拟合C++Templete实现

#include <iostream>
#include <vector>
#include <cmath>
#include <cassert>
#include <cstdlib>

using namespace std;

template<size_t SIZE, size_t N>
class Regression
{
public:
    //ctor
    Regression()
    {
        Y.resize ( SIZE );
        X.resize ( SIZE * N );
        Result.resize ( SIZE );
        SSQ = 0.0L;
    }

    //dtor
    ~Regression() {}
    //get Y[SIZE]
    void _y ( const std :: vector<long double>& y )
    {
        assert ( SIZE == y.size() );
        Y = y;
    }

    //set vector x[][]
    void _x ( const std :: vector<long double>& x )
    {
        assert ( SIZE * N == x.size() );
        X = x;
    }

    //set vector x[][] cloumn index
    void _x ( const std :: vector<long double>& x, const unsigned int& index )
    {
        assert ( index < N );
        assert ( SIZE == x.size() );

        for ( unsigned int i = 0; i < SIZE; ++i )
        {
            X[ i * N + index ] = x[i];
        }
    }

    //execute least square fit
    void fit()
    {
        //the model is
        //             alfa[N][N] * result[N] = beta[N]
        //
        long double alfa[N][N];
        long double beta[N];

        unsigned int i = 0;
		unsigned int j = 0;
		unsigned int k = 0;
		for ( k = 0; k < N; ++k )
        {
            for ( j = 0; j < N; ++j )
            {
                long double tmp = 0.0L;
                for ( i = 0; i < SIZE; ++i )
                {
                    tmp += X[i*N+k] * X[i*N+j];
                }//end of i loop
                alfa[k][j] = tmp;
            }//end of j loop

            long double tmp = 0.0L;
            for ( i = 0; i < SIZE; ++i )
            {
                tmp += Y[i] * X[i*N+k];
            }//end of i loop
            beta[k] = tmp;
        }//end of k loop

        std::vector<long double> a;
        a.clear();
        for ( i = 0; i < N; ++i )
            for ( j = 0; j < N; ++j )
            {
                a.push_back ( alfa[i][j] );
            }

        std::vector<long double> b;
        b.clear();
        for ( i = 0; i < N; ++i )
            b.push_back ( beta[i] );

        //calculate Result[N]
        lq<long double> ( a, b, N, Result );

        //calculate chi-square
        long double chi = 0.0L;
        for ( i = 0; i < SIZE; ++i )
        {
            long double tmp = 0.0L;
            for ( j = 0; j < N; ++j )
            {
                tmp += Result[j] * X[i*j+j];
            }
            tmp -= Y[i];
            chi += tmp * tmp;
        }

        SSQ = gammq ( static_cast<long double> ( SIZE - 2 ) / 2.0L, chi / 2.0L );


    }//end of k loop

    std :: vector< long double> result() const
    {
        return Result;
    }

    long double ssq() const
    {
        return SSQ;
    }

private:
    std :: vector<long double> Y;       // -- Y[SIZE]
    std :: vector<long double> X;       // -- X[N][SIZE]
    std :: vector<long double> Result;  // -- Result[N]
    long double SSQ;                      //fitness of the model
};


static long double
gammln ( const long double &xx )
{
    long double     x,
    tmp,
    ser;
    static long double cof[6] = { 76.18009173, -86.50532033, 24.01409822,
                                  -1.231739516, 0.120858003e-2, -0.536382e-5
                                };
    int             j;

    x = xx - 1.0;
    tmp = x + 5.5;
    tmp -= ( x + 0.5 ) * log ( tmp );
    ser = 1.0;
    for ( j = 0; j <= 5; j++ )
    {
        x += 1.0;
        ser += cof[j] / x;
    }
    return -tmp + log ( 2.50662827465 * ser );
}

static void
gcf ( long double *gammcf,
      const long double &a, const long double &x, long double *gln )
{
    const unsigned int ITMAX = 100;
    const long double EPS = 3.0e-7;
    unsigned int    n = 0;
    long double     gold = 0.0,
                           g,
                           fac = 1.0,
                                 b1 = 1.0;
    long double     b0 = 0.0,
                         anf,
                         ana,
                         an,
                         a1,
                         a0 = 1.0;

    *gln = gammln ( a );
    a1 = x;
    for ( n = 1; n <= ITMAX; n++ )
    {
        an = ( float ) n;
        ana = an - a;
        a0 = ( a1 + a0 * ana ) * fac;
        b0 = ( b1 + b0 * ana ) * fac;
        anf = an * fac;
        a1 = x * a0 + anf * a1;
        b1 = x * b0 + anf * b1;
        if ( a1 )
        {
            fac = 1.0 / a1;
            g = b1 * fac;
            if ( fabs ( ( g - gold ) / g ) < EPS )
            {
                *gammcf = exp ( -x + a * log ( x ) - ( *gln ) ) * g;
                return;
            }
            gold = g;
        }
    }
    if ( "a too large, ITMAX too small in routine GCF" )
        exit ( EXIT_FAILURE );
}

static void
gser ( long double *gamser,
       const long double &a, const long double &x, long double *gln )
{
    const unsigned int ITMAX = 100;
    const long double EPS = 3.0e-7;
    unsigned int    n = 0;
    long double     sum,
    del,
    ap;

    *gln = gammln ( a );
    if ( x <= 0.0 )
    {
        if ( x < 0.0 )
            if ( "x less than 0 in routine GSER" )
                exit ( EXIT_FAILURE );
        *gamser = 0.0;
        return;
    }
    else
    {
        ap = a;
        del = sum = 1.0 / a;
        for ( n = 1; n <= ITMAX; n++ )
        {
            ap += 1.0;
            del *= x / ap;
            sum += del;
            if ( fabs ( del ) < fabs ( sum ) * EPS )
            {
                *gamser = sum * exp ( -x + a * log ( x ) - ( *gln ) );
                return;
            }
        }
        if ( "a too large, ITMAX too small in routine GSER" )
            exit ( EXIT_FAILURE );

        return;
    }
}

long double
gammq ( const long double &a, const long double &x )
{
    long double     gamser,
    gammcf,
    gln;

    if ( x < 0.0 || a <= 0.0 )
        if ( "Invalid arguments in routine gammq." )
            exit ( EXIT_FAILURE );

    if ( x < ( a + 1.0 ) )
    {
        gser ( &gamser, a, x, &gln );
        return 1.0 - gamser;
    }
    else
    {
        gcf ( &gammcf, a, x, &gln );
        return gammcf;
    }
}

//
//for solution of the linear equation
//                  A X = B
//
template<class T>
void lq ( const std :: vector<T>& a,    //store A[size][size]
          const std :: vector<T>& b,     //store B[size]
          const unsigned int order, //store size
          std :: vector<T>& result       //store X[size]
        )
{
    unsigned int Order = order;
    //if order is set to default, calculate
    if ( 0 == Order )
    {
        Order = b.size();
    }

    assert ( Order * Order == a.size() );
    assert ( Order == b.size() );

    result.clear();

    const T _D = det ( a, Order );  //store D

    if ( 0 == _D )
    {
        if ( "Failed to solve the linear equation" )
            exit ( EXIT_FAILURE );
    }

    for ( unsigned int i = 0; i < Order; ++i )
    {
        std :: vector<T> A = a;
        for ( unsigned int j = 0; j < Order; ++j )
        {
            A[i+j*Order] = b[j];
        }
        const T D = det ( A, Order );   //store D[i]

        result.push_back ( D / _D );    //get X[i]
    }
}

//return the determinant of a square matrix arr whose size is order by order
template< class T >                                                         
T det ( const std :: vector<T>& arr, const unsigned int& order = 0 )        
{                                                                           
    unsigned int Order = order;                                             
                                                                            
    //if order is set to default, calculate it                              
    if ( 0 == Order )                                                       
    {                                                                       
        Order = sqrt ( arr.size() );                                        
    }                                                                       
                                                                            
    assert ( Order * Order == arr.size() );                                 
                                                                            
                                                                            
    if ( 1 == Order )                                                       
        return ( arr[0] ) ;                                                 
                                                                            
    T sum = 0;                                                              
    std ::vector<T> arr2 ;                                                  
    int sign = 1;                                                           
                                                                            
    for ( unsigned int i = 0 ; i < Order ; ++i, sign *= -1 )                
    {                                                                       
        /* copy n-1 by n-1 array into another array */                      
        const unsigned int newsize = ( Order - 1 ) * ( Order - 1 ) ;        
                                                                            
        arr2.resize ( newsize );                                            
        for ( unsigned int j = 1 ; j < Order ; ++j )                        
        {                                                                   
            for ( unsigned int k = 0, count = 0 ; k < Order ; ++k )         
            {                                                               
                if ( k == i )                                               
                    continue ;//jump out of k loop                          
                                                                            
                const unsigned int pos = j * Order + k ;                    
                const unsigned int newpos = ( j - 1 ) *                     
                                            ( Order - 1 ) + count ;         
                arr2[newpos] = arr[pos] ;                                   
                count++ ;                                                   
            }//end of k loop                                                
        }//end of j loop                                                    
                                                                            
        /*  find determinant value of n-1 by n-1 array and  add it to sum */
        sum += arr[i] * sign * det ( arr2, Order - 1 ) ;                    
    }                                                                       
    return sum;                                                             
}                                                                           

//下边给出一个示例
//要拟合的参数方程为
//a1 x1 + a2 x2 + a3 x3 + a4 x4 + a5 x5 = y
//给出5组数据进行拟合，分别是（1 0 0 0 0 1）（0 1 0 0 0 2）（0 0 1 0 0 3）（0 0 0 1 0 4）（0 0 0 0 1 5）
int main()
{
    const unsigned int M = 5;
	const unsigned int N = 5;

    Regression<M, N> r;
    vector<long double> y;
    vector<long double> x;
    vector<long double> result;


	long double XXX[N] = { 14.0, 15.7, 17.6, 20.1, 24.7 };
	long double YYY[N] = { 1.541, 1.641, 1.795, 1.668, 1.595 };
    unsigned int i = 0;

	for ( i = 0; i < N; ++i )
        y.push_back( 1.0L + i );
    r._y( y );

    for ( i = 0; i < N; ++i )
    {
        x.clear();
        for ( unsigned int j = 0; j < N; ++j )
        {
            if ( i == j )
                x.push_back( 1.0L );
            else
                x.push_back( 0.0L );
        }
        r._x( x, i );
    }
/*
	for ( i = 0; i < N; ++i )
        y.push_back( YYY[i] );
    r._y( y );

    for ( i = 0; i < N; ++i )
    {
        x.clear();
        for ( unsigned int j = 0; j < N; ++j )
        {
            x.push_back( XXX[j] );
        }
        r._x( x, i );
    }
*/
    r.fit();

    result = r.result();

    for ( i = 0; i < N; ++i )
        cout << result[i] << endl;


    return 0;
}