C++ Math: Fitting a complicated function is not working
I've been trying to get this fitting operation to work since more than 10 days. I wrote a C++ class for fitting using Levenberg-Marquardt algorithm, and I tested it with simple polynomials and the fit is a success; but with the real experimental function I need for my experiment, it's not working.

The same function I fitted with Mathematica, and it was OK to do it there (but slow, 12 seconds per fit, which is why I'm using C++). My program acts crazy and returns wrong results when I do the fit. Could you please check the functions and derivatives I use and tell me whether anything is wrong?

The following is the function I wanna fit and its parameters' derivatives. SR is the sampling rate, it's a given constant; and Real type is a typedef for double:


    Real FIDFunc(Real t, Real Intercept, Real AmpExp,  Real TauExp, Real AmpMod, Real TauMod, Real Freq, Real Phase)
    {
        return Intercept + AmpExp*exp(-t/SR/TauExp) + AmpMod*exp(-t/SR/TauMod)*cos((2*M_PI*Freq*t)/SR + Phase);
    }
    Real FIDFunc_D_Intercept(Real t, Real Intercept, Real AmpExp,  Real TauExp, Real AmpMod, Real TauMod, Real Freq, Real Phase)
    {
        return 1;
    }
    Real FIDFunc_D_AmpExp(Real t, Real Intercept, Real AmpExp,  Real TauExp, Real AmpMod, Real TauMod, Real Freq, Real Phase)
    {
        return exp(-t/SR/TauExp);
    }
    Real FIDFunc_D_TauExp(Real t, Real Intercept, Real AmpExp,  Real TauExp, Real AmpMod, Real TauMod, Real Freq, Real Phase)
    {
        return AmpExp*t/(SR*TauExp*TauExp)*exp(-t/SR/TauExp);
    }
    Real FIDFunc_D_AmpMod(Real t, Real Intercept, Real AmpExp,  Real TauExp, Real AmpMod, Real TauMod, Real Freq, Real Phase)
    {
        return exp(-t/SR/TauMod)*cos((2*M_PI*Freq*t)/SR + Phase);
    }
    Real FIDFunc_D_TauMod(Real t, Real Intercept, Real AmpExp,  Real TauExp, Real AmpMod, Real TauMod, Real Freq, Real Phase)
    {
        return AmpMod*t/(SR*TauMod*TauMod)*exp(-t/SR/TauMod)*cos((2*M_PI*Freq*t)/SR + Phase);
    }
    Real FIDFunc_D_Freq(Real t, Real Intercept, Real AmpExp,  Real TauExp, Real AmpMod, Real TauMod, Real Freq, Real Phase)
    {
        return -AmpMod*2*M_PI*t/SR*exp(-t/SR/TauMod)*sin((2*M_PI*Freq*t)/SR + Phase);
    }
    Real FIDFunc_D_Phase(Real t, Real Intercept, Real AmpExp,  Real TauExp, Real AmpMod, Real TauMod, Real Freq, Real Phase)
    {
        return -AmpMod*exp(-t/SR/TauMod)*sin((2*M_PI*Freq*t)/SR + Phase);
    }

The following is the procedure of fitting which I use in my class. Again, I'd like to mention that it worked with polynomials, so I could've have done a very deep mistake here that conflicts with only some special functions... that's my assumption.


    template <typename T>
    const SamArray<T>& SamFit<T>::fit()
    {
        T lambda = 0.001;
        SamArray<T> parameters = params; //params is a member variable of the class, which is an array of the parameters
        T prevChiSquare;
    
        for(unsigned long i = 0; i < MAXITERATIONS; i++)
        {
            getChiSquare(params);
            getJacobian(params);
            prevChiSquare = ChiSquare; //ChiSquare is sum of the squares of the differences. It's a member variable in the class too
            parameters = params;
            SamArray<SamArray<T> > combinationMatrix = SamArray<T>::DiagonalMatrix(JTJ); //JTJ is a member variable, which is the Kronecker product of the Jacobian with itself.
            for(long j = 0; j < combinationMatrix.size(); j++)
            {
                combinationMatrix[j] *= lambda; //here I'm following the Wikipedia gradient acceleration
                //check http://en.wikipedia.org/wiki/Levenberg%E2%80%93Marquardt_algorithm
            }
            parameters += SamArray<T>::MatrixProduct(SamArray<T>::Inverse(JTJ + combinationMatrix),J); //add the new deltas to the parameters and test them
    
            getChiSquare(parameters);
            if(ChiSquare < prevChiSquare)
            {
                params = parameters; //accept the updated parameters
                lambda /= 10;
            }
            else
            {
                parameters = params; //reject the updated parameters
                lambda *= 10;
            }
            if(abs(ChiSquare - prevChiSquare) < tolerance)
            {
                cout<<"Break..."<<endl;
                break;
            }
        }
        return params;
    }

The initial values I use are very close to the correct ones.

If you need any additional information, please let me know. Any help is very highly appreciated.

Thank you.
C++ Math: Fitting a complicated function is not working
===

I've been trying to get this fitting operation to work since more than 10 days. I wrote a C++ class for fitting using Levenberg-Marquardt algorithm, and I tested it with simple polynomials and the fit is a success; but with the real experimental function I need for my experiment, it's not working.

The same function I fitted with Mathematica, and it was OK to do it there (but slow, 12 seconds per fit, which is why I'm using C++). My program acts crazy and returns wrong results when I do the fit. Could you please check the functions and derivatives I use and tell me whether anything is wrong?

The following is the function I wanna fit and its parameters' derivatives. SR is the sampling rate, it's a given constant; and Real type is a typedef for double:


    Real FIDFunc(Real t, Real Intercept, Real AmpExp,  Real TauExp, Real AmpMod, Real TauMod, Real Freq, Real Phase)
    {
        return Intercept + AmpExp*exp(-t/SR/TauExp) + AmpMod*exp(-t/SR/TauMod)*cos((2*M_PI*Freq*t)/SR + Phase);
    }
    Real FIDFunc_D_Intercept(Real t, Real Intercept, Real AmpExp,  Real TauExp, Real AmpMod, Real TauMod, Real Freq, Real Phase)
    {
        return 1;
    }
    Real FIDFunc_D_AmpExp(Real t, Real Intercept, Real AmpExp,  Real TauExp, Real AmpMod, Real TauMod, Real Freq, Real Phase)
    {
        return exp(-t/SR/TauExp);
    }
    Real FIDFunc_D_TauExp(Real t, Real Intercept, Real AmpExp,  Real TauExp, Real AmpMod, Real TauMod, Real Freq, Real Phase)
    {
        return AmpExp*t/(SR*TauExp*TauExp)*exp(-t/SR/TauExp);
    }
    Real FIDFunc_D_AmpMod(Real t, Real Intercept, Real AmpExp,  Real TauExp, Real AmpMod, Real TauMod, Real Freq, Real Phase)
    {
        return exp(-t/SR/TauMod)*cos((2*M_PI*Freq*t)/SR + Phase);
    }
    Real FIDFunc_D_TauMod(Real t, Real Intercept, Real AmpExp,  Real TauExp, Real AmpMod, Real TauMod, Real Freq, Real Phase)
    {
        return AmpMod*t/(SR*TauMod*TauMod)*exp(-t/SR/TauMod)*cos((2*M_PI*Freq*t)/SR + Phase);
    }
    Real FIDFunc_D_Freq(Real t, Real Intercept, Real AmpExp,  Real TauExp, Real AmpMod, Real TauMod, Real Freq, Real Phase)
    {
        return -AmpMod*2*M_PI*t/SR*exp(-t/SR/TauMod)*sin((2*M_PI*Freq*t)/SR + Phase);
    }
    Real FIDFunc_D_Phase(Real t, Real Intercept, Real AmpExp,  Real TauExp, Real AmpMod, Real TauMod, Real Freq, Real Phase)
    {
        return -AmpMod*exp(-t/SR/TauMod)*sin((2*M_PI*Freq*t)/SR + Phase);
    }

The following is the procedure of fitting which I use in my class. Again, I'd like to mention that it worked with polynomials, so I could've have done a very deep mistake here that conflicts with only some special functions... that's my assumption.


    template <typename T>
    const SamArray<T>& SamFit<T>::fit()
    {
        T lambda = 0.001;
        SamArray<T> parameters = params; //params is a member variable of the class, which is an array of the parameters
        T prevChiSquare;
    
        for(unsigned long i = 0; i < MAXITERATIONS; i++)
        {
            getChiSquare(params);
            getJacobian(params);
            prevChiSquare = ChiSquare; //ChiSquare is sum of the squares of the differences. It's a member variable in the class too
            parameters = params;
            SamArray<SamArray<T> > combinationMatrix = SamArray<T>::DiagonalMatrix(JTJ); //JTJ is a member variable, which is the Kronecker product of the Jacobian with itself.
            for(long j = 0; j < combinationMatrix.size(); j++)
            {
                combinationMatrix[j] *= lambda; //here I'm following the Wikipedia gradient acceleration
                //check http://en.wikipedia.org/wiki/Levenberg%E2%80%93Marquardt_algorithm
            }
            parameters += SamArray<T>::MatrixProduct(SamArray<T>::Inverse(JTJ + combinationMatrix),J); //add the new deltas to the parameters and test them
    
            getChiSquare(parameters);
            if(ChiSquare < prevChiSquare)
            {
                params = parameters; //accept the updated parameters
                lambda /= 10;
            }
            else
            {
                parameters = params; //reject the updated parameters
                lambda *= 10;
            }
            if(abs(ChiSquare - prevChiSquare) < tolerance)
            {
                cout<<"Break..."<<endl;
                break;
            }
        }
        return params;
    }

The initial values I use are very close to the correct ones.

If you need any additional information, please let me know. Any help is very highly appreciated.

Thank you.