Matlab options for curve-fitting with custom equation (TDR application for RC measurement)

Question

I have the following custom function that must be fitted to the data:

function y=reflectometriaRC(x,v0,t0,v1,t1,v2,v3,tau)

y=zeros(size(x));

for i=1:length(x)
    if x(i)<t0,
        y(i)=v0;
    elseif x(i)<t0+t1,
        y(i)=v1;
    else
        y(i)=v2+v3*(1-exp(-(x(i)-(t0+t1))/tau));
    end
end
end

You can see an example of this function in the following image:

The problem is that i cannot find the correct options to make the "fit" function to work propertly.

I have tried things like this:

x=transpose([0:200*10^-12:1200*10^-9]);
y=reflectometriaRC(x,0,100*10^-9,0.98,100*10^-9,0,1.96,250*10^-9);
ft = fittype('reflectometriaRC(x,v0,t0,v1,t1,v2,v3,tau)');
f = fit( x, y, ft,'Robust','on', 'MaxFunEvals',600000,'MaxIter',400000000, 'StartPoint', [0, 0, 0.5, 0, 0,0, 0], 'Lower',[0, 0, 0.5, 0, 0,0, 0]  );
coeffvalues(f)

But the result is just some sort of step function.

Any idea what would be a good configuration to make the fitting work?

---

The idea of this code is to get a measure from an oscilloscope used with a pulse generator like a basic TDR. The device under test in an RC circuit at the end of a 1m coaxial line. The image in the link above is in the scale of the possible values in seconds and volt.

Answer 1

In general, fitting a multi-parameter model to underconstrained data is a hard problem. You might ultimately have to use more sophisticated approaches, such as using a Markov chain Monte Carlo approach -- which has the benefit of being guaranteed to converge to the global error minimum (albeit not in a predictable amount of time).

Have you tried using a different fitting algorithm? If you've got access to the statistics toolbox, try doc nlinfit. This function can use either the Levenberg-Marquardt algoirthm or an iterative reweighted least squares algorithm. Moreover, you can get confidence intervals relatively easily for your fitted parameters, which gives you a good idea about how good the fit is.

The reason fit doesn't work very well is that it, by default, uses a Levenberg-Marquardt algorithm which relies upon gradient descent methods that are prone to getting stuck in local minima. So, to be fair, is nlinfit, but you've got a bit more control over it.

First of all, redefine your function slightly:

function y=reflectometriaRC(b, x)
%The function to fit; inputs: b, a vector of parameters, and x, 
%the independent variable (e.g. time)
v0=b(1);
t0=b(2);
v1=b(3);
t1=b(4);
v2=b(5);
v3=b(6);
tau=b(7);

y=zeros(size(x));

for i=1:length(x)
    if x(i)<t0,
        y(i)=v0;
    elseif x(i)<t0+t1,
        y(i)=v1;
    else
        y(i)=v2+v3*(1-exp(-(x(i)-(t0+t1))/tau));
    end
end
end

Then, in your script:

opts = statset('nlinfit');
opts.RobustWgtFun = 'bisquare';
beta0= [0, 0, 0.5, 0, 0,0, 0];
params=nlinfit(x, y, @reflectometricaRC, beta0,opts);

This returns parameters that are relatively close to the true ones -- numerically, three of the four are identical (zero) and the others are within a factor of two, but the fit still isn't great. You're likely going to either have to play around with beta0, the error model used, or -- my advice -- use a lot more computational time and go for the MCMC approach -- it's what I use when fitting high-dimensional models.

Hope that helps!