Fitting a sigmoid curve using curve fitting tool box

Hi, I am trying to fit a sigmoid function to the underlying data with the goodness of fit. I was using the curve fitting tool box. Unfortunately, i am not getting an idea to how to go forward with this. I tried to write custom function to fit the data, unfortunately, i am not able to get the required fit. Any help to solve this will be appreciated. X and Y data is shown below.

y= [0.70,1.09,0.96,0.37,0.27,0.35,0.16,0.18,0.18,0.02,0.37,0.67,1.26,1.70,1.85,1.73,1.92,1.72,1.47,1.50,1.46,1.43,2.12,1.6,1.8]

x = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25]

Answer 1

[Alex Sha]

Hi, Hair, the quailty of your data seems to be bad for sigmoid type of function. If take the function like:

y = y0+a/(1+exp(-(x-x0)/b))^c

Then the fitting result looks like:

[Hari krishnan]

@ Alexshah.

Thank you very much. Can you explain to me why did you choose this equation and what does each term in the equation implies?

[Alex Sha]

Hi, Hari, since you did not given the clear expression of fitting function, but only emphasized on "Sigmoid", so I just selected a Sigmoid-type function, there are many such kind of function avaiable.

[Hari krishnan]

@ Alex Sha, the data corresponds to population dynamics. So the data captures the trend of a regular logistic growth. I wanted to know what the value of coefficient signifies?

Answer 2

[John D'Errico]

The problem is, there are many functional forms that people think of and call with the generic name "sigmoid". We can find a few of them listed here:

https://en.wikipedia.org/wiki/Sigmoid_function

Perhaps the most common one is of the specific form:

f(x) = 1/(1 + exp(-x))

But many people will want to use thing like a transformed cumulative normal, which can be gotten from the erf function. Or possibly an arc tangent curve, or many others. Again, a few of them are shown on the wiki page. The thing is, each of those has a different fundamental shape.

y = [0.70,1.09,0.96,0.37,0.27,0.35,0.16,0.18,0.18,0.02,0.37,0.67,1.26,1.70,1.85,1.73,1.92,1.72,1.47,1.50,1.46,1.43,2.12,1.6,1.8]
x = 1:length(y);
plot(x,y,'o')

Your first problem is that your data is REALLY noisy. The first three points almost look like they are not even from the same curve.

Regardless, lets try to fit a simple sigmoidal model to your data. I'll start with the one I show above, adding some parameters, but keeping the same general form.

ft = fittype('L + (U-L)/(1 + exp(-4*log(3)*(x-xmid)/xscale80))','indep','x')
ft = 
     General model:
     ft(L,U,xmid,xscale80,x) = L + (U-L)/(1 + exp(-4*log(3)*(x-xmid)/xscale80))
     
mdl = fit(x',y',ft,'start',[.5 ,3, 15,3])
mdl = 
     General model:
     mdl(x) = L + (U-L)/(1 + exp(-4*log(3)*(x-xmid)/xscale80))
     Coefficients (with 95% confidence bounds):
       L =      0.4242  (0.2428, 0.6057)
       U =       1.694  (1.519, 1.869)
       xmid =       12.69  (11.91, 13.47)
       xscale80 =       1.763  (-1.042, 4.567)

Here, I've named the parameters to be indicative of what they represent in the model. As it turns out, each of the paramters can be simply interpreted here.

That is, L is the lower asymptote, the limit of f(x), as x --> inf. U is the upper asymptote. xmid is the point where the curve goes through 50% of the transition. So xmid is the 50% point.

xscale80 is a little more difficult, but we can express it in the form of the width of the transition. If my arithmetic is correct, we can say the difference between the 10% and the 90% points on the curve will be xscale80. So 80% of the transition occurs between the 10% nd 90% points on that curve. Thus, when xscale80 is a large number, the transition is a slow one. If xscale is small, then the transition approaches a sharp step function.

plot(mdl)
grid on
hold on
plot(x,y,'bo')
xline(mdl.xmid);
xline(mdl.xmid + mdl.xscale80/2);
xline(mdl.xmid - mdl.xscale80/2);

I've drawn vertical lines at the 10% 50% and 90% points on the curve as it goes through the transition.

As you can see, the model does fit quasi-reasonably, although the lower asymptote is a bit off in my eyes. I did say your data was pure crapola, I think? At least those first three data points are dragging the estimate of L off from where it seems it should be. A weighted fit, or a robust regression would help.

mdl = fit(x',y',ft,'Robust','LAR','start',[.5 ,3, 15,3])
Fitting stopped because the number of iterations or function evaluations exceeded the specified maximum.
mdl = 
     General model:
     mdl(x) = L + (U-L)/(1 + exp(-4*log(3)*(x-xmid)/xscale80))
     Coefficients (with 95% confidence bounds):
       L =      0.2993  (0.1037, 0.4948)
       U =        1.72  (1.532, 1.908)
       xmid =       12.59  (11.74, 13.44)
       xscale80 =       2.374  (-0.7873, 5.536)
       
plot(mdl)
hold on
grid on
plot(x,y,'bo')
xline(mdl.xmid - mdl.xscale80/2);
xline(mdl.xmid);
xline(mdl.xmid + mdl.xscale80/2);

This fit looks more reasonable. Your data is still pure crapola though. This is hinted at by the width of the confidence intervals on the parameters. As you can see, those bounds on the parameters are all scarily wide.

[Hari krishnan]

thanks!

[Srilatha Raghavan]

Hey John, I have a quick question in continutaion of this discussion. I am using a sigmoid curve to fit my power curve(Wind speed Vs Power). So I need to force the curve to stay at 0 and start only at, 3 m/s on the X-axis and abrubtly stop at 25m/s on X-axis. I can never seem to custom write a code for this because then my sigmoid is behaving like off. Could you help me identify the component that can be adjusted for this?? Thank you in advance.

In short, I want to force the fitted curve to stay at 0 for a while before fitting the points. However, my data need not be 0 during that period.

[John D'Errico]

If you want to fit a sigmoid to data, you need to understand that the curve you have chosen does not switch on and off like a light bulb. Any such sigmoidal model will only appoach the lower and upper asymptotes, but these are asymptotes.

x = sort(rand(1,50)*40)
y = 3*(erf((x - 22)/5)+1) + randn(size(x))/5;
plot(x,y,'o')
xon = 10;
xoff = 32;
xline(xon);
xline(xoff);

I've constructed a curve that has the indicated basic shape, but then I'll pretend I need to fit it. Beelow the first line, I'll force the function to be identically 0. Above that point, I'll force it to be constant, though I won't specify any known value.

There are two simple ways you might choose to do such a fit. Either use my SLM toolbox, as found on the file exchange for download. Or you could use a truncated sigmoidal function.

slm = slmengine(x,y,'plot','on','knots',[0 10 13:4:29 32 40],'increasing','on','leftvalue',0,'leftslope',0,'xyp',[10,0;32,0],'rightslope',0,'concaveup',[0,18],'concavedown',[25,40]);
xline(xon);
xline(xoff);

The fit seems quite good. It has all the properties you would want to see. It seems to estimate the upper asymptote well enough.

slmeval(40,slm)
ans =
          6.05814929210623

Alternatively, I could use the curve fitting toolbox, with a custom model. I'll use an m-file to create it.

function yhat = mysigmoid(a,b,c,x,xon,xoff)
  x = max(xon,min(x,xoff));
  yhat = a./(1 + exp(-(x - c)/b)) - a./(1 + exp(-(xon - c)/b));
end

Now call fit.

xon = 10;
xoff = 32;
ft = fittype(@(a,b,c,x) mysigmoid(a,b,c,x,xon,xoff));

mdl = fit(x',y',ft,'start',[5  5 20])
mdl = 
     General model:
     mdl(x) = mysigmoid(a,b,c,x,xon,xoff)
     Coefficients (with 95% confidence bounds):
       a =       6.173  (6.016, 6.329)
       b =       2.224  (1.974, 2.474)
       c =       22.08  (21.85, 22.32)
       
figure
plot(mdl)
hold on
plot(x,y,'bo')
xline(xon);
xline(xoff);

Again, a very reasonable fit. Below x == 10, it is identically zero. About 32 it is constant.

mdl([0;10;32;40])
ans =
                         0
                         0
          6.07511093779466
          6.07511093779466

[Srilatha Raghavan]

Hello. Thank you so much for getting back to me. Your code works perfectly fine. I still have one more thing that is taking a long time to figure out if possible.

I want my fitted curve to extend till, say x=25. I wrote a five line code which pulls my fitted curve from 0 in x-axis and stops at x=25.

In short, I want my fitted curve to start from 0 and go on till 25. Is that possible?

[John D'Errico]

I'm not sure what you are asking. The code I wrote allows the function to be evaluated for any value of x. But why is that a problem? If you don't want to evaluate it above or below those limits, then don't try to do so.

Are you asking the function to produce an error if you try to pass it a value of x below 0 or above 25? Personally, that seems a bit silly. Writing code that will intentionally result in an error when you do something reasonable seems wrong. Instead, The code I wrote returns the maximum or minimum value above or below the supplied limits.

[Srilatha Raghavan]

Sorry I just realized. That was a stupid question. I can just high pass filter my input data.

Just one last thing. Is it possible to fix the value of a. Because when I do so, its not being adjusted. When I try to remove the variable 'a' out of the picture and have my number in the equation(I don't wanna do this), I get this error:

Custom equations must produce an output vector, matrix, or array that is the same size and shape as the input data. This custom equation fails to meet that requirement.

I wanna fix 'a' just like 'xon' and 'xoff'. This is basically something down to a very stupid mistake I am doing.

 xon = 3;
 xrated = 12.5;
 xoff = 25;
 a =900;
 ft = fittype(@(b,c,x) mysigmoid(b,c,x,xon,xrated,a));
 mdl = fit(x',y',ft,'start',[5 20]);
 
 function yhat = mysigmoid(b,c,xon,xrated,a)
 x= max(xon,min(x,xrated));
 yhat = a./(1+exp(-(x-c)/b)) - a./(1+exp(-(xon-c)/b));
 end

TIA

[Srilatha Raghavan]

It's the sigmoid function ( the equation) adjustment I am having troubles about. There's always a trade off between the sigmoid curve I got from your equation and my equation adjustment to get my curve to reach a particular height, at my rated x value.