头文件:
/*
* 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 <iostream>
#include <constants.h>
namespace splab
{
template <typename Type>
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 <matrix.h>
#include <linequs1.h>
#include <fitcurves.h>
#include <interpolation.h>
namespace splab
{
template <typename Type>
class LSFitting : public Interpolation<Type>
{
public:
using Interpolation<Type>::xi;
using Interpolation<Type>::yi;
LSFitting( const Vector<Type> &xn, const Vector<Type> &yn,
Funcs<Type> &f );
~LSFitting();
void calcCoefs();
Type evaluate( Type x );
Vector<Type> getCoefs() const;
private:
Vector<Type> coefs; // fitted parameters
Funcs<Type> phi; // fitted functions
};
// class LSFitting
#include <lsfitting-impl.h>
}
// 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 <typename Type>
LSFitting<Type>::LSFitting( const Vector<Type> &xn, const Vector<Type> &yn,
Funcs<Type> &f )
: Interpolation<Type>( xn, yn ), phi(f)
{
coefs.resize(phi.dim);
}
template <typename Type>
LSFitting<Type>::~LSFitting()
{
}
/**
* Compute fitted parameters.
*/
template <typename Type>
void LSFitting<Type>::calcCoefs()
{
int N = xi.size(),
M = coefs.dim();
Type tmp;
Vector<Type> b(M);
Matrix<Type> C( M, M );
for( int i=1; i<=M; ++i )
{
tmp = 0;
for( int k=0; k<N; ++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<N; ++k )
tmp += phi(i,xi[k]) * phi(j,xi[k]);
C(i,j) = C(j,i) = tmp;
}
tmp = 0;
for( int k=0; k<N; ++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 <typename Type>
Type LSFitting<Type>::evaluate( Type x )
{
Type y = 0;
for( int j=0; j<coefs.size(); ++j )
y += coefs[j] * phi( j, x );
return y;
}
/**
* Get the fitted parameters.
*/
template <typename Type>
inline Vector<Type> LSFitting<Type>::getCoefs() const
{
return coefs;
}
测试代码:
/*****************************************************************************
* lsfit_test.cpp
*
* Least square fitting testing.
*
* Zhang Ming, 2010-04, Xi'an Jiaotong University.
*****************************************************************************/
#define BOUNDS_CHECK
#include <iostream>
#include <lsfitting.h>
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<Type> x = linspace( 0.01, 0.5*PI, M );
Vector<Type> y(M);
for( int i=0; i<M; ++i )
{
tmp = A * pow(x[i],alpha) * exp(beta*x[i]+gamma*x[i]*x[i]);
y[i] = log(max(tmp,EPS));
}
Funcs<Type> phi;
LSFitting<Type> lsf( x, y, phi );
lsf.calcCoefs();
Vector<Type> 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.