Why is the fit() function producing a horizontal line above my data?

I am using the MATLAB fit() function (single exponential) to fit the decay of miniature excitaotry post synaptic current (mEPSC) events. It works really nicely for some events, but for others the fitted curve is simply a flat line above my data (see represemtative image for an example of each). I cannot empirically determine any differences between the events where the fit looks good compared to those where the fit produces the horizontal line. Any ideas?

[Walter Roberson]

Do I understand correctly that you are extracting x and y ranges from your data, and fitting each of the segments with exp1 model, ?

(That model can be computed easily with the \ operator after taking log of both sides.)

[Luke Fournier]

Yes, the starting point is the peak (the minimum), the end is when it returns to baseline. This is all determined programatically in my code. And yes, I am using exp1 with the equation you listed. Thank you for your prompt reply.

Any ideas?

[Walter Roberson]

Could you extract the data for that segment near 37 and attach the x and y coordinates here, so that we could test ?

[Luke Fournier]

Here is an excel that has the data. The first two columns (as labeled) have a larger segment of the trace around 37, the right most two columns have the x and y data for the exact segment I am attemtping to fit. Thank you.

[Torsten]

I hope you don't work with x-values of about 35000 in the fitting function a*exp(b*x) ...

[Luke Fournier]

This is a segment of a longer trace and it happens to start at data point ~35000. The current is sampled every 0.05ms for 60s, so naturally the whole trace has quite a large number of data points. However, what the actual data I am trying to fit consists of only about 130 data points. In this way, the x value of 35000 corresponds to about time point 1.75s out of the 60s trace. For this reason, I am confused as to why the value of the x axis matters.

Perhaps you can clarify. Thank you!

[Luke Fournier]

To clarify further, I could have simply extracted these points and made the first point I gave you 1 and the final data point would have an x value of about 135. In the original picture I attached, you can see the data fit nicely, and that starting x value for the section to fir corresponds to roughly x=30000. I hope this helps

[Luke Fournier]

37, for example, indicates that this event is the 37th that my code identified in entire trace

[the cyclist]

@Torsten, large values of x are OK there, because they can be scaled by the values of b. (But one does need to be careful to specify a good starting value for b to avoid the blow-up that I assume you were cautioning against.)

[the cyclist]

@Luke Fournier, is your model really of the form a*exp(b*x), as @Walter Roberson asked?

Unless I got myself confused, that model cannot fit the data you posted. exp(b*x) is positive everywhere, so your fitted response values will either all be positive, or all negative (depending on the sign of a). Because your response variable has both positive and negative values, you cannot get a good fit.

[Luke Fournier]

@the cyclist thank you for your input. I believe that is the correct form of the equation, as I am certain that I used fit(x,y, ‘exp1’), which utilizes that equation (at least that is my understanding). This is why I am confused as to why it can fit some events but other times it yields the flat line.

[Luke Fournier]

So I know that model works for some of the events, so this is my confusion as to why with some events which look similar in nature it fails.

[Luke Fournier]

If I provided the data for event 36 (as I did for event 37), would this help? So you can see the data fitting for one event but not the other? The code I have written is what identifies which points the take to fit (i.e. the blue data points), but maybe if I provided the data of the raw trace of an event where the fit did work (36) would help? Let me know!

Answer 1

[Walter Roberson]

Activate the "center and scale" option -- but remember that the coefficients you are shown as a result refer to the centered/scaled coordinate system.

[Luke Fournier]

I will try this and update once I do, thanks Walter!

[Luke Fournier]

Problem completely resolved, thank you Walter! For those who may read in the future and wonder about the exact syntax that corresponds to Walter's solution, use fit(x, y, 'exp1', 'Normalize', 'on'). Activating the normalization is what centers/scales the coordinate system. Read more here about fitoptions

Thank you all again for your help!

[Walter Roberson]

Using @Torsten's flow template from https://www.mathworks.com/matlabcentral/answers/424895-how-can-i-make-logaritmic-fitting-like-polyfit-function#comment_624760 :

format long g
filename = 'for_walter.xlsx';
T = readmatrix(filename);

fx = T(:,4);
fy = T(:,5);
%data contains lots of nan
mask = ~isfinite(fx) | ~isfinite(fy);
fx(mask) = [];
fy(mask) = [];

p = polyfit(fx,log(fy),1)

p = 
         -0.023959908827361 -    0.0167790652408582i           847.540790585289 +      595.387462252473i

b = p(1);
a = p(2);

syms x
yfit = exp(sym(a)) * exp(b*x);
for_display = vpa((yfit),16)

for_display = 

Notice the coefficients on the order of 6e+366 -- coefficients that are beyond the range representable by double precision. So when you work in double precision, you cannot get a useful fit.

Let us try a different approach:

syms areal aimag breal bimag
a = areal + 1i*aimag;
b = breal + 1i*bimag;
%residue = sum((fy - a*exp(b*fx)).^2);
residue = norm(fy - a*exp(b*fx));
rfun = matlabFunction(residue, 'vars', {[areal, aimag, breal, bimag]});
[best, fval] = fmincon(rfun, [-.1 -.1 1e-5 1e-5], [], [], [], [], [-1e6 -1e6 -700 -700], [1e6 1e6 700 700])

Local minimum possible. Constraints satisfied.

fmincon stopped because the size of the current step is less than
the value of the step size tolerance and constraints are 
satisfied to within the value of the constraint tolerance.

best = 1×4

1.0e+00 * -359574.445229496 -148362.537586897 -0.000324026729102402 -0.000189093661138137

fval = 
          45.7651903184081

a = best(1) + 1i*best(2)

a = 
          -359574.445229496 -      148362.537586897i

b = best(3) + 1i*best(4)

b = 
      -0.000324026729102402 -  0.000189093661138137i

y = a*exp(b*x)

y = 

vpa(y)

ans = 

scatter(fx, fy);
hold on
fplot(real(y), [min(fx), max(fx)])
hold off

Coefficients not quite as large, but the plot is still flat.

If you try to fit a*exp(b*x) using fmincon and using only real values, Yes, you get an answer, but the residue is large