forum.alglib.net
http://forum.alglib.net/

levenberg-marquardt lsfit with and without weights
http://forum.alglib.net/viewtopic.php?f=2&t=1405
Page 1 of 1

Author:  asadullahbeg [ Sat Feb 01, 2014 8:15 am ]
Post subject:  levenberg-marquardt lsfit with and without weights

Hi,

I am implementing a function as seen in the attachment. the parameters which are changing are a,b and c0. I have the x and y data from experimentation but sadly I am not getting any correct values of a,b and c. I am still not implementing the example with the hessian and gradient included. Should I implement that example? Is there anything wrong with my code? Please help :(
Attachment:
equation.png
equation.png [ 3.19 KiB | Viewed 9621 times ]








Code:
#include "stdafx.h"
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#define PI 3.1415926535897
#include "interpolation.h"

using namespace alglib;
void function_cx_1_func(const real_1d_array &c, const real_1d_array &x, double &func, void *ptr)
{
    // this callback calculates f(c,x)=exp(-c0*sqr(x0))
    // where x is a position on X-axis and c is adjustable parameter
   func = (1 / (c[0] * c[1])) * (c[2] - c[0] * 1 * sinh(sinh(c[2]) * exp(c[0] * x[0] * 2 *PI)));
        //func= 1/ab*(c-a*sinh(sinh(c))*exp(a*x*2*pi);
    //func = exp(-c[0]*pow(x[0],2));
}
void function_cx_1_grad(const real_1d_array &c, const real_1d_array &x, double &func, real_1d_array &grad, void *ptr)
{
    // this callback calculates f(c,x)=exp(-c0*sqr(x0)) and gradient G={df/dc[i]}
    // where x is a position on X-axis and c is adjustable parameter.
    // IMPORTANT: gradient is calculated with respect to C, not to X
   
   double temp;
     grad[0] = -(c[2]);
            grad[0] = grad[0] / (pow(c[0], 2) * c[1]);
            temp=(2*PI*x[0]*exp(2*PI*c[0]*x[0])*cosh(exp(2*PI*c[0]*x[0])*sinh(c[2]))*sinh(c[2]));
           temp=temp/c[1];
            grad[0]=grad[0]-temp;
           grad[1] = sinh( exp(2 * PI * c[0] * x[0]) *  sinh(c[2])) /  pow(c[1], 2) - c[2] / (c[0] *  pow(c[1], 2));
            grad[2] = 1 / (c[0] * c[1]) - ( exp(2 *  PI * c[0] * x[0]) *  cosh( exp(2 *  PI * c[0] * x[0]) *  sinh(c[2])) *  cosh(c[2])) / c[1];

}
void function_cx_1_hess(const real_1d_array &c, const real_1d_array &x, double &func, real_1d_array &grad, real_2d_array &hess, void *ptr)
{
    // this callback calculates f(c,x)=exp(-c0*sqr(x0)), grad G={df/dc[i]} and Hessian H={d2f/(dc[i]*dc[j])}
    // where x is a position on X-axis and c is adjustable parameter.
    // IMPORTANT: gradient/Hessian are calculated with respect to C, not to X
    func = exp(-c[0]*pow(x[0],2));
   hess[0][0]=(2*c[2])/(pow(c[0],3)*c[1])-(4*pow(PI,2)*pow(x[0],2)*exp(4*PI*c[0]*x[0])*sinh(exp(2*PI*c[0]*x[0])*sinh(c[2]))*pow(sinh(c[2]),2))/c[1] - (4*pow(PI,2)*pow(x[0],2)*exp(2*PI*c[0]*x[0])*cosh(exp(2*PI*c[0]*x[0])*sinh(c[2]))*sinh(c[2]))/c[1];
   hess[0][1]=c[2]/(pow(c[0],2)*pow(c[1],2)) + (2*PI*x[0]*exp(2*PI*c[0]*x[0])*cosh(exp(2*PI*c[0]*x[0])*sinh(c[2]))*sinh(c[2]))/pow(c[1],2);
   hess[0][2]= - 1/(pow(c[0],2)*c[1]) - (2*PI*x[0]*exp(2*PI*c[0]*x[0])*cosh(exp(2*PI*c[0]*x[0])*sinh(c[2]))*cosh(c[2]))/c[1] - (2*PI*x[0]*exp(4*PI*c[0]*x[0])*sinh(exp(2*PI*c[0]*x[0])*sinh(c[2]))*cosh(c[2])*sinh(c[2]))/c[1];
   hess[1][0]= c[2]/(pow(c[0],2)*pow(c[1],2)) + (2*PI*x[0]*exp(2*PI*c[0]*x[0])*cosh(exp(2*PI*c[0]*x[0])*sinh(c[2]))*sinh(c[2]))/pow(c[1],2);
   hess[1][1]= (2*c[2])/(c[0]*pow(c[1],3)) - (2*sinh(exp(2*PI*c[0]*x[0])*sinh(c[2])))/pow(c[1],3);
   hess[1][2]=(exp(2*PI*c[0]*x[0])*cosh(exp(2*PI*c[0]*x[0])*sinh(c[2]))*cosh(c[2]))/pow(c[1],2) - 1/(c[0]*pow(c[1],2));
   hess[2][0]= - 1/(pow(c[0],2)*c[1]) - (2*PI*x[0]*exp(2*PI*c[0]*x[0])*cosh(exp(2*PI*c[0]*x[0])*sinh(c[2]))*cosh(c[2]))/c[1] - (2*PI*x[0]*exp(4*PI*c[0]*x[0])*sinh(exp(2*PI*c[0]*x[0])*sinh(c[2]))*cosh(c[2])*sinh(c[2]))/c[1];
   hess[2][1]= (exp(2*PI*c[0]*x[0])*cosh(exp(2*PI*c[0]*x[0])*sinh(c[2]))*cosh(c[2]))/pow(c[1],2) - 1/(c[0]*pow(c[1],2));
   hess[2][2]= - (exp(2*PI*c[0]*x[0])*cosh(exp(2*PI*c[0]*x[0])*sinh(c[2]))*sinh(c[2]))/c[1] - (exp(4*PI*c[0]*x[0])*sinh(exp(2*PI*c[0]*x[0])*sinh(c[2]))*pow(cosh(c[2]),2))/c[1];
   
}

int main(int argc, char **argv)
{
    //
    // In this example we demonstrate exponential fitting
    // by f(x) = exp(-c*x^2)
    // using function value, gradient and Hessian (with respect to c)
    //
    real_2d_array x = "[[1],[2],[3],[4],[5],[6],[7]]";
   
    real_1d_array y = "[438,469,563,797,1196,1410,1721]";
    real_1d_array c = "[0.0,0.22,0.0]";
     double epsf = 0;
    double epsx = 0.000001;
    ae_int_t maxits = 0;
    ae_int_t info;
    lsfitstate state;
    lsfitreport rep;
    double diffstep = 0.0001;

    //
    // Fitting without weights
    //
    lsfitcreatef(x, y, c, diffstep, state);
    lsfitsetcond(state, epsf, epsx, maxits);
    alglib::lsfitfit(state, function_cx_1_func);
    lsfitresults(state, info, c, rep);
    printf("%d\n", int(info)); // EXPECTED: 2
    printf("%s\n", c.tostring(1).c_str()); // EXPECTED: [1.5]

    //
    // Fitting with weights
    // (you can change weights and see how it changes result)
    //
    real_1d_array w = "[1,1,1,1,1,1,1]";
    lsfitcreatewf(x, y, w, c, diffstep, state);
    lsfitsetcond(state, epsf, epsx, maxits);
    alglib::lsfitfit(state, function_cx_1_func);
    lsfitresults(state, info, c, rep);
   
    printf("%s\n", c.tostring(1).c_str()); // EXPECTED: [1.5](v
   (void)getchar();
    return 0;
}

Author:  asadullahbeg [ Sun Feb 02, 2014 5:16 am ]
Post subject:  Re: levenberg-marquardt lsfit with and without weights

I have implemented this routine for other functions and experimental data all comes out smooth. Maybe my experimental data is wrong? Maybe the estimated function is wrong? Please I need help :(

Author:  Sergey.Bochkanov [ Sun Feb 02, 2014 7:42 am ]
Post subject:  Re: levenberg-marquardt lsfit with and without weights

Hello! I see following problems:

1. you start at point A=0, but your function involves division by A. Division by zero produces infinity which confuses optimization algorithm. You should start from another point, non-zero A and B. It is better to guard A and B away from zero by means of boundary constraints, like A>0.00001, B>0.00001.

2. it is better to set non-zero initial c0.

Author:  asadullahbeg [ Tue Feb 04, 2014 3:50 am ]
Post subject:  Re: levenberg-marquardt lsfit with and without weights

real_1d_array c = "[0.0525,0.075,-1.35]";

Thank you Sergey that solved the problem. Ofc I had to find these values with many trials. I had a question if the initial guess isn't right meaning real_1d_array c = "[1,1,-1]"; will the LM algo still get the right a,b and c0? Because I wasnt getting the right values of the a,b and c0. How I know that? Because when i placed these values in the function and saw the plot the answer came to be a straight line. Again thanks a million.

Author:  Sergey.Bochkanov [ Tue Feb 04, 2014 7:11 am ]
Post subject:  Re: levenberg-marquardt lsfit with and without weights

If I understand you right, it is highly problem specific. For some problems, LM can converge from almost any point. For other ones, it needs good initial point. Your problem seems to be hard one - complex shape, easy to stuck into local extrema.

If you want something like 100% automated solution:
1) set boundary constraints on A and B, as I recommended you
2) run LM algorithm several times from random A/B/C0, and choose best result

Author:  asadullahbeg [ Sat Feb 08, 2014 9:57 am ]
Post subject:  Re: levenberg-marquardt lsfit with and without weights

Hi Sergey,
New problem. How can I get the values of real_1d_array from user?

Author:  asadullahbeg [ Sat Feb 08, 2014 10:00 am ]
Post subject:  Re: levenberg-marquardt lsfit with and without weights

Hi Sergey,
New problem. How can I get the values of real_1d_array from user?

Author:  asadullahbeg [ Sat Feb 08, 2014 10:06 am ]
Post subject:  Re: levenberg-marquardt lsfit with and without weights

using c++

Author:  Sergey.Bochkanov [ Sun Feb 09, 2014 9:31 am ]
Post subject:  Re: levenberg-marquardt lsfit with and without weights

"get the values of real_1d_array from user" - what does it mean? :)

Page 1 of 1 All times are UTC
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group
http://www.phpbb.com/