最小二乘曲线拟合 C语言算法实现一,最小二乘曲线拟合算法的C++实现

头文件：

/*

* Copyright (c) 2008-2011 Zhang Ming (M. Zhang), zmjerry@163.com

*

* This program is free software; you can redistribute it and/or modify it

* under the terms of the GNU General Public License as published by the

* Free Software Foundation, either version 2 or any later version.

*

* Redistribution and use in source and binary forms, with or without

* modification, are permitted provided that the following conditions are met:

*

* 1. Redistributions of source code must retain the above copyright notice,

* this list of conditions and the following disclaimer.

*

* 2. Redistributions in binary form must reproduce the above copyright

* notice, this list of conditions and the following disclaimer in the

* documentation and/or other materials provided with the distribution.

*

* This program is distributed in the hope that it will be useful, but WITHOUT

* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or

* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for

* more details. A copy of the GNU General Public License is available at:

* http://www.fsf.org/licensing/licenses

*/

/*****************************************************************************

* fitcurves.h

*

* Fitted functions for least square fitting.

*

* Zhang Ming, 2010-04, Xi'an Jiaotong University.

*****************************************************************************/

#ifndef FITCURVES_H

#define FITCURVES_H

#include

#include

namespace splab

{

template

class Funcs

{

public:

const static int dim = 4;

// Compute the value of functions at point x.

Type operator()( int i, const Type &x )

{

switch( i )

{

case 1:

{

return 1;

break;

}

case 2:

{

return log(max(x,EPS));

break;

}

case 3:

{

return x;

break;

}

case 4:

{

return x*x;

break;

}

default:

{

std::cerr << "The dimension 'i' exceed the bound!"

<< std::endl;

return 0;

break;

}

}

}

};

// class ObjFunc

}

// namespace splab

#endif

// FITCURVES_H

/*

* Copyright (c) 2008-2011 Zhang Ming (M. Zhang), zmjerry@163.com

*

* This program is free software; you can redistribute it and/or modify it

* under the terms of the GNU General Public License as published by the

* Free Software Foundation, either version 2 or any later version.

*

* Redistribution and use in source and binary forms, with or without

* modification, are permitted provided that the following conditions are met:

*

* 1. Redistributions of source code must retain the above copyright notice,

* this list of conditions and the following disclaimer.

*

* 2. Redistributions in binary form must reproduce the above copyright

* notice, this list of conditions and the following disclaimer in the

* documentation and/or other materials provided with the distribution.

*

* This program is distributed in the hope that it will be useful, but WITHOUT

* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or

* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for

* more details. A copy of the GNU General Public License is available at:

* http://www.fsf.org/licensing/licenses

*/

/*****************************************************************************

* lsfitting.h

*

* Least Square Fitting.

*

* For a given points set "Pn" {xn,yn} and a fitted function with paramter

* "pm", this class can find the best parameter "pm" in the meaning of least

* mean-square error.

*

* Zhang Ming, 2010-04, Xi'an Jiaotong University.

*****************************************************************************/

#ifndef LSFITTING_H

#define LSFITTING_H

#include

#include

#include

#include

namespace splab

{

template

class LSFitting : public Interpolation

{

public:

using Interpolation::xi;

using Interpolation::yi;

LSFitting( const Vector &xn, const Vector &yn,

Funcs &f );

~LSFitting();

void calcCoefs();

Type evaluate( Type x );

Vector getCoefs() const;

private:

Vector coefs; // fitted parameters

Funcs phi; // fitted functions

};

// class LSFitting

#include

}

// namespace splab

#endif

// LSFITTING_H

实现文件：

/*

* Copyright (c) 2008-2011 Zhang Ming (M. Zhang), zmjerry@163.com

*

* This program is free software; you can redistribute it and/or modify it

* under the terms of the GNU General Public License as published by the

* Free Software Foundation, either version 2 or any later version.

*

* Redistribution and use in source and binary forms, with or without

* modification, are permitted provided that the following conditions are met:

*

* 1. Redistributions of source code must retain the above copyright notice,

* this list of conditions and the following disclaimer.

*

* 2. Redistributions in binary form must reproduce the above copyright

* notice, this list of conditions and the following disclaimer in the

* documentation and/or other materials provided with the distribution.

*

* This program is distributed in the hope that it will be useful, but WITHOUT

* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or

* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for

* more details. A copy of the GNU General Public License is available at:

* http://www.fsf.org/licensing/licenses

*/

/*****************************************************************************

* lsfitting-impl.h

*

* Implementation for LSFitting class.

*

* Zhang Ming, 2010-04, Xi'an Jiaotong University.

*****************************************************************************/

/**

* constructors and destructor

*/

template

LSFitting::LSFitting( const Vector &xn, const Vector &yn,

Funcs &f )

: Interpolation( xn, yn ), phi(f)

{

coefs.resize(phi.dim);

}

template

LSFitting::~LSFitting()

{

}

/**

* Compute fitted parameters.

*/

template

void LSFitting::calcCoefs()

{

int N = xi.size(),

M = coefs.dim();

Type tmp;

Vector b(M);

Matrix C( M, M );

for( int i=1; i<=M; ++i )

{

tmp = 0;

for( int k=0; k

tmp += phi(i,xi[k]) * phi(i,xi[k]);

C(i,i) = tmp;

for( int j=i+1; j<=M; ++j )

{

tmp = 0;

for( int k=0; k

tmp += phi(i,xi[k]) * phi(j,xi[k]);

C(i,j) = C(j,i) = tmp;

}

tmp = 0;

for( int k=0; k

tmp += phi(i,xi[k]) * yi[k];

b(i) = tmp;

}

coefs = choleskySolver( C, b );

}

/**

* Compute the value of fitted function at given "x".

*/

template

Type LSFitting::evaluate( Type x )

{

Type y = 0;

for( int j=0; j

y += coefs[j] * phi( j, x );

return y;

}

/**

* Get the fitted parameters.

*/

template

inline Vector LSFitting::getCoefs() const

{

return coefs;

}

测试代码：

/*****************************************************************************

* lsfit_test.cpp

*

* Least square fitting testing.

*

* Zhang Ming, 2010-04, Xi'an Jiaotong University.

*****************************************************************************/

#define BOUNDS_CHECK

#include

#include

using namespace std;

using namespace splab;

typedef double Type;

int main()

{

Type A = 4.0,

alpha = 1.0,

beta = -2.0,

gamma = -4.0,

tmp = 0.0;

int M = 100;

Vector x = linspace( 0.01, 0.5*PI, M );

Vector y(M);

for( int i=0; i

{

tmp = A * pow(x[i],alpha) * exp(beta*x[i]+gamma*x[i]*x[i]);

y[i] = log(max(tmp,EPS));

}

Funcs phi;

LSFitting lsf( x, y, phi );

lsf.calcCoefs();

Vector parms = lsf.getCoefs();

parms(1) = exp(parms(1));

cout << "The original parameters are:" << endl

<< A << endl << alpha << endl << beta << endl << gamma << endl << endl;

cout << "The fitted parameters are: " << parms << endl;

return 0;

}

运行结果：

The original parameters are:

4

1

-2

-4

The fitted parameters are: size: 4 by 1

4

1

-2

-4

Process returned 0 (0x0) execution time : 0.094 s

Press any key to continue.