필터 지우기
필터 지우기

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

조회 수: 7 (최근 30일)
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')

채택된 답변

John D'Errico
John D'Errico 2024년 4월 6일
편집: John D'Errico 2024년 4월 7일
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.
  댓글 수: 3
John D'Errico
John D'Errico 2024년 4월 7일
편집: John D'Errico 2024년 4월 7일
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
Mushtaq Al-Jubbori 2024년 4월 7일
@John D'Errico I want model with few parameters,excruciatingly well and pass through 0,1 and used for long x=0 to 80

댓글을 달려면 로그인하십시오.

추가 답변 (1개)

John D'Errico
John D'Errico 2024년 4월 6일
편집: John D'Errico 2024년 4월 6일
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')
  댓글 수: 4
Mushtaq Al-Jubbori
Mushtaq Al-Jubbori 2024년 4월 6일
편집: Mushtaq Al-Jubbori 2024년 4월 6일
@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
John D'Errico 2024년 4월 6일
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.

댓글을 달려면 로그인하십시오.

카테고리

Help CenterFile Exchange에서 Linear and Nonlinear Regression에 대해 자세히 알아보기

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by