can any one help me to find simple empirical model (equation) with a little parameters for this data

x=[31.96 30.25 28.27 26.84 25.4 24 23.27 21.89 20.85 18.6 17.02 16 15.81 13.6 12 11.19 9.8 8.8 7.07 5.7 4.6 4.36 3.2 2.6 2.15 1.5 1 0.75 0.35 0.1 ]; y=[1.18 1.26 1.38 1.5 1.63 1.75 1.89 2.068 2.21 2.55 2.8 2.99 3 3.43 3.7 3.89 4.1 4.3 4.7 5 5.2 5.3 5.5 5.6 5.7 5.7 5.4 4.9 3.2 1.16]; plot(x,y,'o')

[Mushtaq Al-Jubbori]

@Alex Sha

Answer 1

[John D'Errico]

As an alternate model, one that employs the assumption that the curve passes through the point (0,1), can we find a different fit? For example, we could try a sum of exponentials, but with that assumption built in.

x=[31.96 30.25 28.27 26.84 25.4 24  23.27 21.89 20.85 18.6 17.02 16 15.81 13.6 12 11.19 9.8 8.8 7.07 5.7 4.6 4.36 3.2 2.6 2.15 1.5 1 0.75 0.35 0.1  ];
y=[1.18 1.26 1.38 1.5 1.63 1.75 1.89 2.068 2.21 2.55 2.8 2.99 3 3.43 3.7 3.89 4.1 4.3 4.7 5 5.2 5.3 5.5 5.6 5.7 5.7 5.4 4.9 3.2 1.16];

mdl = fittype('a*exp(-b*x) + (1-a)*exp(-c*x)',indep = 'x')

mdl = 
     General model:
     mdl(a,b,c,x) = a*exp(-b*x) + (1-a)*exp(-c*x)

fittedmdl = fit(x',y',mdl,start = [6 2 1])

fittedmdl = 
     General model:
     fittedmdl(x) = a*exp(-b*x) + (1-a)*exp(-c*x)
     Coefficients (with 95% confidence bounds):
       a =       6.756  (6.533, 6.979)
       b =     0.05321  (0.05025, 0.05617)
       c =       1.555  (1.364, 1.746)

plot(fittedmdl,x,y,'bo')

And again, we can see the curve now passes exactly through (0,1).

fittedmdl(0)

ans = 1

As well though, you can see the fit is rather poor. I could throw another exponential term in there of course, with the result that now the model gets more complex.

mdl = fittype('a*exp(-b*x) + (c)*exp(-d*x) +  (1-a-c)*exp(-e*x)',indep = 'x')

mdl = 
     General model:
     mdl(a,b,c,d,e,x) = a*exp(-b*x) + (c)*exp(-d*x) +  (1-a-c)*exp(-e*x)

fittedmdl = fit(x',y',mdl,start = [5 2 1 .1 .2])

fittedmdl = 
     General model:
     fittedmdl(x) = a*exp(-b*x) + (c)*exp(-d*x) +  (1-a-c)*exp(-e*x)
     Coefficients (with 95% confidence bounds):
       a =       958.8  (-1.564e+08, 1.564e+08)
       b =       1.167  (-427.1, 429.5)
       c =       5.371  (4.47, 6.272)
       d =     0.03882  (0.02841, 0.04922)
       e =       1.172  (-427, 429.4)

plot(fittedmdl,x,y,'bo')

With total crap for a result, based on wild guesses for starting values.

Is there a valid reason to chase down this rabbit hole?

Ok. For kicks, to make it easy to understand how this was done, I'll use GA, coupled with an optimproblem to estimate the parameters.

prob = optimproblem;
linv = optimvar('linv',[1,3],LowerBound = -10,UpperBound = 10);
ratev = optimvar('ratev',[1,3],LowerBound = 0);
Prob.constraints.MyConstraint1 = sum(linv) == 1;
% the objective function is simple. 
prob.Objective = sum((linv(1)*exp(-ratev(1)*x) + linv(2)*exp(-ratev(2)*x) + linv(3)*exp(-ratev(3)*x) - y).^2);

sol = solve(prob,Solver = 'ga')

Solving problem using ga.
ga stopped because it exceeded options.MaxGenerations.

sol = struct with fields:
     linv: [-3.9585 -6.1918 9.9375]
    ratev: [0.1471 2.4366 0.0680]

plot(x,y,'bo')
hold on
fplot(@(x) exp(-x(:).*sol.ratev)*sol.linv',[0,max(x)])

Warning: Function behaves unexpectedly on array inputs. To improve performance, properly vectorize your function to return an output with the same size and shape as the input arguments.

hold off
grid on

Honestly, that took me a few tries to get even GA to find a decent solution.

[Mushtaq Al-Jubbori]

@John D'Errico May be coplex model but stisfay the condutions at x=0, y=1

[John D'Errico]

The other models I also show, ALSO satisy that condition! And all of them fit your data far better. Only the very first one I built did not use that condition, BECAUSE you never told me that was a requirement.

Make up your mind! Do you want a fit with few parameters? Do you want a fit that fits excruciatingly well, with many parameters? It seems that you are asking for both of these things, at different times.

I understand now that you want a fit that passes through the point (0,1), but I have shown many ways to do that fit. Do you want a fit that can be extrapolated for large x? Or do you only care what it does within the support of your data?

[Mushtaq Al-Jubbori]

@John D'Errico I want model with few parameters,excruciatingly well and pass through 0,1 and used for long x=0 to 80

Answer 2

[John D'Errico]

Why not try it? This seems adequate:

x=[31.96 30.25 28.27 26.84 25.4 24  23.27 21.89 20.85 18.6 17.02 16 15.81 13.6 12 11.19 9.8 8.8 7.07 5.7 4.6 4.36 3.2 2.6 2.15 1.5 1 0.75 0.35 0.1  ];
y=[1.18 1.26 1.38 1.5 1.63 1.75 1.89 2.068 2.21 2.55 2.8 2.99 3 3.43 3.7 3.89 4.1 4.3 4.7 5 5.2 5.3 5.5 5.6 5.7 5.7 5.4 4.9 3.2 1.16];

So I'll put a negative exponential in there for the fast transient in the beginning, then a simple polynomial that will take over when x grows large.

mdl = fittype('a + b*x + c*x^2 + d*x.^3+ e*exp(-f*x)',indep = 'x')

mdl = 
     General model:
     mdl(a,b,c,d,e,f,x) = a + b*x + c*x^2 + d*x.^3+ e*exp(-f*x)

fittedmdl = fit(x',y',mdl,start = [1 1 .1 .1 -6 2])

fittedmdl = 
     General model:
     fittedmdl(x) = a + b*x + c*x^2 + d*x.^3+ e*exp(-f*x)
     Coefficients (with 95% confidence bounds):
       a =       6.221  (6.165, 6.277)
       b =     -0.2184  (-0.2334, -0.2033)
       c =  -3.266e-06  (-0.001047, 0.00104)
       d =   5.919e-05  (3.85e-05, 7.988e-05)
       e =       -6.35  (-6.446, -6.253)
       f =       2.259  (2.183, 2.335)

plot(fittedmdl,x,y,'bo')

[Mushtaq Al-Jubbori]

Dear @John D'Errico Thank you very much. But if you can found another equation with 3 or 4 parameters. You can add (1) to the equation to reduced the parameters

[John D'Errico]

@Mushtaq Al-Jubbori

Sigh. Ok. So do it. But sadly, that won't reduce the number of parameters! Think about what you see, before jumping to conclusions.

Note that the constant coefficient as estimated was 6.221, NOT 1. The difference seems significant to me.

And, yes, you COULD have said that you had information about the model? That the curve should pass exactly through the point (0,1)? Double sigh. I hate chasing a moving target.

If you do know that, AND you want to enforce it, then the model should change. Not by simply adding 1. What happens in the model I used when x==0, and I try to force the curve through that point? We would get this:

a + e == 1

Change the model to reflect that.

x=[31.96 30.25 28.27 26.84 25.4 24  23.27 21.89 20.85 18.6 17.02 16 15.81 13.6 12 11.19 9.8 8.8 7.07 5.7 4.6 4.36 3.2 2.6 2.15 1.5 1 0.75 0.35 0.1  ];
y=[1.18 1.26 1.38 1.5 1.63 1.75 1.89 2.068 2.21 2.55 2.8 2.99 3 3.43 3.7 3.89 4.1 4.3 4.7 5 5.2 5.3 5.5 5.6 5.7 5.7 5.4 4.9 3.2 1.16];
mdl = fittype('a + b*x + c*x^2 + d*x.^3 + (1-a)*exp(-f*x)',indep = 'x')

mdl = 
     General model:
     mdl(a,b,c,d,f,x) = a + b*x + c*x^2 + d*x.^3 + (1-a)*exp(-f*x)

fittedmdl = fit(x',y',mdl,start = [6 1 .1 .1 2])

fittedmdl = 
     General model:
     fittedmdl(x) = a + b*x + c*x^2 + d*x.^3 + (1-a)*exp(-f*x)
     Coefficients (with 95% confidence bounds):
       a =       6.401  (6.065, 6.737)
       b =     -0.2564  (-0.3425, -0.1703)
       c =    0.002228  (-0.003635, 0.008091)
       d =   2.032e-05  (-9.455e-05, 0.0001352)
       f =       1.714  (1.463, 1.965)

plot(fittedmdl,x,y,'bo')

fittedmdl(0)

ans = 1

So now it passes EXACTLY through the point (0,1). And that was implicitly what you asked to see. However, do you now see the fit is considerably worse? I could probably try using a different model of course, one that has extra parameters. But that would seem palpably silly, since you wanted to REDUCE the number of parameters.

In the original model, as I fit it, you would see that the curve fits well, but it does not really want to pass through the point (0,1). If that seems important to you, then you need to use a more complex model, not a simpler one. So, is that now the goal? Are there any other factors I need to take into considereration? For example, do you want to examine the long term behavior of this empirically derived model? Remember, you asked only for an empirical model that fits the data!

[Mushtaq Al-Jubbori]

@John D'Errico Thanks, when I asked, I say with a little parameters. The equation with allitle parameters was beeter Nessary pass through (0,1) When x=0, must y=1

[John D'Errico]

Sigh. Exactly how do the words "little parameters" imply the curve must pass through (0,1)?

You say that you want FEWER parameters. I want to see world peace in my lifetime. Does my wanting something make it come true? Why should that be true for you?

We can add extra terms of course.

x=[31.96 30.25 28.27 26.84 25.4 24  23.27 21.89 20.85 18.6 17.02 16 15.81 13.6 12 11.19 9.8 8.8 7.07 5.7 4.6 4.36 3.2 2.6 2.15 1.5 1 0.75 0.35 0.1  ];
y=[1.18 1.26 1.38 1.5 1.63 1.75 1.89 2.068 2.21 2.55 2.8 2.99 3 3.43 3.7 3.89 4.1 4.3 4.7 5 5.2 5.3 5.5 5.6 5.7 5.7 5.4 4.9 3.2 1.16];
mdl = fittype('a + b*x + c*x^2 + d*x.^3 + e*x.^4 + (1-a)*exp(-f*x)',indep = 'x')

mdl = 
     General model:
     mdl(a,b,c,d,e,f,x) = a + b*x + c*x^2 + d*x.^3 + e*x.^4 + (1-a)*exp(-f*x)

fittedmdl = fit(x',y',mdl,start = [6 -0.2 0.002 0.0001 0.000001 1.5])

fittedmdl = 
     General model:
     fittedmdl(x) = a + b*x + c*x^2 + d*x.^3 + e*x.^4 + (1-a)*exp(-f*x)
     Coefficients (with 95% confidence bounds):
       a =       6.638  (6.042, 7.234)
       b =     -0.3599  (-0.5822, -0.1376)
       c =     0.01476  (-0.0102, 0.03972)
       d =  -0.0005378  (-0.001611, 0.0005351)
       e =   8.205e-06  (-7.35e-06, 2.376e-05)
       f =        1.63  (1.342, 1.919)

plot(fittedmdl,x,y,'bo')

That seems to follow the data better in some places, but it also may be chasing noise in the curve too.

I wold go with the model in this form, as you saw in my last comment.

mdl = fittype('a + b*x + c*x^2 + d*x.^3 + (1-a)*exp(-f*x)',indep = 'x')

It satisfies the requirement that is passes through (0,1), but it is not chasing noise in your data. Higher order models do that.

If you want a simpler model, for gods sake, you asked for an empirical model. Simple does not always equate to good in terms of fit. In fact, very often you will find a tradeoff. Simpler models are aften less good models.