pythonmatlaboptimizationcurve-fittinglevenberg-marquardt

Reducing difference between two graphs by optimizing more than one variable in MATLAB/Python?


Suppose 'h' is a function of x,y,z and t and it gives us a graph line (t,h) (simulated). At the same time we also have observed graph (observed values of h against t). How can I reduce the difference between observed (t,h) and simulated (t,h) graph by optimizing values of x,y and z? I want to change the simulated graph so that it imitates closer and closer to the observed graph in MATLAB/Python. In literature I have read that people have done same thing by Lavenberg-marquardt algorithm but don't know how to do it?


Solution

  • You are actually trying to fit the parameters x,y,z of the parametrized function h(x,y,z;t).

    MATLAB

    You're right that in MATLAB you should either use lsqcurvefit of the Optimization toolbox, or fit of the Curve Fitting Toolbox (I prefer the latter).

    Looking at the documentation of lsqcurvefit:

    x = lsqcurvefit(fun,x0,xdata,ydata);
    

    It says in the documentation that you have a model F(x,xdata) with coefficients x and sample points xdata, and a set of measured values ydata. The function returns the least-squares parameter set x, with which your function is closest to the measured values.

    Fitting algorithms usually need starting points, some implementations can choose randomly, in case of lsqcurvefit this is what x0 is for. If you have

    h = @(x,y,z,t) ... %// actual function here
    t_meas = ... %// actual measured times here
    h_meas = ... %// actual measured data here
    

    then in the conventions of lsqcurvefit,

    fun   <--> @(params,t) h(params(1),params(2),params(3),t)
    x0    <--> starting guess for [x,y,z]: [x0,y0,z0]
    xdata <--> t_meas
    ydata <--> h_meas
    

    Your function h(x,y,z,t) should be vectorized in t, such that for vector input in t the return value is the same size as t. Then the call to lsqcurvefit will give you the optimal set of parameters:

    x = lsqcurvefit(@(params,t) h(params(1),params(2),params(3),t),[x0,y0,z0],t_meas,h_meas);
    h_fit = h(x(1),x(2),x(3),t_meas);  %// best guess from curve fitting
    

    Python

    In python, you'd have to use the scipy.optimize module, and something like scipy.optimize.curve_fit in particular. With the above conventions you need something along the lines of this:

    import scipy.optimize as opt
    
    popt,pcov = opt.curve_fit(lambda t,x,y,z: h(x,y,z,t), t_meas, y_meas, p0=[x0,y0,z0])
    

    Note that the p0 starting array is optional, but all parameters will be set to 1 if it's missing. The result you need is the popt array, containing the optimal values for [x,y,z]:

    x,y,z = popt
    h_fit = h(x,y,z,t_meas)