I need to fits the attached data as in image
Show older comments
I need to fit the attached data as in the image below:
2 Comments
Scott MacKenzie
on 5 Jun 2021
Have you made any attempt yourself to do this? If no, then perhaps you should. If yes, then it might help if you showed the code you've put together so far.
Mushtaq Al-Jubbori
on 5 Jun 2021
Accepted Answer
Alex Sha
on 6 Jun 2021
Hi, the fitting function "y = p1/(1+Exp(p2+p3*x)+p4*x)^p5" is pretty good for all data set, where y=L, x=T
1: T1-L1:
Root of Mean Square Error (RMSE): 0.0481024350848356
Sum of Squared Residual: 0.0277661311330899
Correlation Coef. (R): 0.999085841980414
R-Square: 0.998172519645713
Parameter Best Estimate
---------- -------------
p1 6.26675447680147
p2 109.105635295479
p3 -67.1792841156261
p4 404935160.81852
p5 0.0161845071634169
2: T2-L2:
Root of Mean Square Error (RMSE): 0.0728301819569392
Sum of Squared Residual: 0.111388943481498
Correlation Coef. (R): 0.999334306206347
R-Square: 0.998669055560921
Parameter Best Estimate
---------- -------------
p1 7.42506477062961
p2 19.1556613661477
p3 -9.26234046573527
p4 -0.0219401101448161
p5 0.101288835010354
3: T3-L3:
Root of Mean Square Error (RMSE): 0.134092463526771
Sum of Squared Residual: 0.413558141817605
Correlation Coef. (R): 0.999276494987953
R-Square: 0.998553513435409
Parameter Best Estimate
---------- -------------
p1 11.4516493147292
p2 17.1095166703901
p3 -4.66159981573592
p4 0.0192821696183541
p5 0.129232672490336
4: T4-L4:
Root of Mean Square Error (RMSE): 0.174646688455142
Sum of Squared Residual: 1.06755130259216
Correlation Coef. (R): 0.999289186413636
R-Square: 0.998578878083226
Parameter Best Estimate
---------- -------------
p1 15.0727185425809
p2 61.2119358895463
p3 -10.7637556382665
p4 0.236072827848106
p5 0.037477078210382
5: T5-L5:
Root of Mean Square Error (RMSE): 0.154312409495508
Sum of Squared Residual: 0.500058714210495
Correlation Coef. (R): 0.999654251549009
R-Square: 0.999308622640009
Parameter Best Estimate
---------- -------------
p1 17.5989211228526
p2 138.340503536049
p3 -17.2818282854456
p4 -0.0357446437398955
p5 0.0179266207350469
11 Comments
Image Analyst
on 6 Jun 2021
Looks like that model fits pretty well. Where does that model come from? Is it the theoretical model for some actual real world physical process?
Scott MacKenzie
on 7 Jun 2021
@Alex Sha Could you please post the code you used to generate one of the above fits and plots.
Alex Sha
on 7 Jun 2021
Hi, the function was got automatically through another package of 1stOpt. Actually, the original function of "y=A1*tanh(U)^A4, where U=exp((x-A2)/A3)" is also good enough:
1: T1-L1:
Root of Mean Square Error (RMSE): 0.0597205155644429
Sum of Squared Residual: 0.0427984797513945
Correlation Coef. (R): 0.998591934579049
R-Square: 0.997185851806327
Parameter Best Estimate
---------- -------------
a1 4.49354483105404
a2 1.2952779431858
a3 0.182566357732861
a4 0.204961211278745
2: T2-L2:
Root of Mean Square Error (RMSE): 0.0787759184709166
Sum of Squared Residual: 0.130318551949667
Correlation Coef. (R): 0.999222163002948
R-Square: 0.998444931036289
Parameter Best Estimate
---------- -------------
a1 7.50386733656648
a2 2.02357393772506
a3 0.425853364237124
a4 0.421653184517995
3: T3-L3:
Root of Mean Square Error (RMSE): 0.130459406471193
Sum of Squared Residual: 0.391452104946764
Correlation Coef. (R): 0.99931463290777
R-Square: 0.99862973554359
Parameter Best Estimate
---------- -------------
a1 11.2849810802541
a2 3.62185746086777
a3 0.583978847866495
a4 0.356432637089129
4: T4-L4:
Root of Mean Square Error (RMSE): 0.17944782733443
Sum of Squared Residual: 1.12705329572666
Correlation Coef. (R): 0.999239734136879
R-Square: 0.998480046277941
Parameter Best Estimate
---------- -------------
a1 14.4986426026826
a2 5.50830739927991
a3 0.561356452592117
a4 0.235099297830832
5: T5-L5:
Root of Mean Square Error (RMSE): 0.185048265050244
Sum of Squared Residual: 0.719100068360212
Correlation Coef. (R): 0.999502208306672
R-Square: 0.999004664409913
Parameter Best Estimate
---------- -------------
a1 17.7699852315831
a2 8.00653745200654
a3 0.42767787817109
a4 0.133811896038511
Scott MacKenzie
on 7 Jun 2021
Edited: Scott MacKenzie
on 7 Jun 2021
Am I missing something? I don't see any code in your comment.
Alex Sha
on 7 Jun 2021
The below is the code for T5-L5:
Variable x,y;
Function y=A1*tanh(exp((x-A2)/A3))^A4;
Data;
x=[2:0.5:7.5 7.7 8:0.5:10.5 11.5 12];
y=[3 3.4 3.8 4.4 5.2 6 6.8 8 9.2 10.8 12.8 15 16.4 17.4 17.8 17.6 17.8 17.8 17.6 17.8 17.8];
when running in 1stOpt, the result will be:
Root of Mean Square Error (RMSE): 0.185048265050244
Sum of Squared Residual: 0.719100068360216
Correlation Coef. (R): 0.999502208325959
R-Square: 0.999004664448469
Parameter Best Estimate
---------- -------------
a1 17.7700127061749
a2 8.00656197680192
a3 0.427679189766609
a4 0.133811270063284
Image Analyst
on 7 Jun 2021
OK, but we don't have 1stOpt - I don't even know what it is. Is there some code that runs in MATLAB? With standard fitting functions like exist in base MATLAB or the Stats Toolbox?
Sulaymon Eshkabilov
on 8 Jun 2021
Edited: Image Analyst
on 9 Jun 2021
It seems that 1stOpt is a 3rd party package or application. It would be nice if @Alex Sha could elaborate more on 1stOpt that demonstrates to be quite accurate. Thanks.
Mushtaq Al-Jubbori
on 10 Jun 2021
Image Analyst
on 10 Jun 2021
Well, I gave you code for that here:
You're free to develop the function further if you want.
Mushtaq Al-Jubbori
on 10 Jun 2021
More Answers (5)
Sulaymon Eshkabilov
on 6 Jun 2021
0 votes
You can start using curve fitting toolbox, cftool that is quite straightforward and does not require any addional coding.
An alternative way is the least squares method for that you have you data ready. You'd need to generate Vandermonde matrix and then compute the fit model coefficients.
4 Comments
Image Analyst
on 6 Jun 2021
What is the Vandermonde model? Do you have any formula for it? If so we might be able to use the Statistics and Machine Learning Toolbox, which I believe is more common than the Curve Fitting Toolbox. Or you could pick some other of various sigmoid functions, such as the Michaelis-Menton rate equation. It would be good to know the theory behind what generated the data so we could pick a proper theoretical equation.
Sulaymon Eshkabilov
on 6 Jun 2021
Vandermonde matrix can be built based on many different formulations. Its beauty and advantage is that it is computationally very efficient. E.g: Fit model: ax^2+bx+c
% y is [Measured DATA];
x = linspace(0, 5); x =x';
V = [x.^2 x ones(size(x))];
Fit_Model = V\y
John D'Errico
on 6 Jun 2021
NO. You do NOT want to use a polynomial model for that data!
Polynomials NEVER have the property that they roll over and then become flat.
Sulaymon Eshkabilov
on 8 Jun 2021
@John D'Errico This polynomial model is shown here mean to be an example not the proposed model. For V matrix it is viable to insert exp(), sin(), cos() tan() and polynonials in combination.
Image Analyst
on 6 Jun 2021
@Mushtaq Al-Jubbori. What I did was to find where the flat part starts and fit the data to the left of that to an exponential growth curve. As an output you have the location where the flat part starts and the coefficients of the exponential fitted equation. Try this (click the copy icon in the upper right of the code block below and then paste into MATLAB):
clc; % Clear command window.
fprintf('Running %s.m ...\n', mfilename);
clear; % Delete all variables.
close all; % Close all figure windows except those created by imtool.
workspace; % Make sure the workspace panel is showing.
fontSize = 14;
L1=[1.4 1.9 2.4 3.132 4.2 4.532 4.532 4.4 4.532 4.532 4.4 4.5];%1.5Mev
L2=[1.65 2.2 2.8 3.4 4.266 5.6 6.6 7.4 7.5 7.534 7.5 7.5 7.4 7.4 7.4 7.5 7.6 7.532 7.4 7.6 7.6];%2.26MeV
L3=[1.8 2.2 2.8 3 3.8 4.8 5.8 6.4 7.8 8.6 9.8 10.6 11.2 11.4 11.2 11.4 11.2 11.4 11 11.2 11.4 11.2 11.4 ];%3.2MeV
L4=[2.05 2.4 2.8 3.4 3.7 4 4.5 5 6.3 7.2 8 8.6 9.8 10.4 11.2 12.2 13.8 14.6 14.6 14.532 14.532 14.6 14.6 14.532 14.532 14.4 14.4 14.4 14.4 14.532 14.5 14.2 14.4 14.4 14.6];%4MeV
L5=[3 3.4 3.8 4.4 5.2 6 6.8 8 9.2 10.8 12.8 15 16.4 17.4 17.8 17.6 17.8 17.8 17.6 17.8 17.8];%4.85MeV
T1=0.25:0.25:3;
T2=[0.5:0.25:2.25 2.75:0.25:5 5.5 5.75 6];
T3=[0.75:0.25:1.75 2.25:0.25:4 4.5:0.25:5 5.5 6 6.25 6.75:0.25:7.5 ];
T4=[1 1.25 1.75:0.25:3 3.5:0.25:9.25 9.75 10.25 10.75];
T5=[2:0.5:7.5 7.7 8:0.5:10.5 11.5 12];
subplot(5, 2, 1);
plot(T1,L1)
grid on;
title('L1 vs. T1', 'FontSize', fontSize);
kneeIndex = FindKneeIndex(L1, T1);
xline(T1(kneeIndex), 'Color', 'r', 'LineWidth', 2);
% Now fit the left portion of the data that follows an exponential to an exponential formula.
coefficients1 = FitExponentialGrowth(T1(1:kneeIndex), L1(1:kneeIndex), 2)
subplot(5, 2, 3);
plot(T2,L2)
grid on;
title('L2 vs. T2', 'FontSize', fontSize);
kneeIndex = FindKneeIndex(L2, T2);
xline(T2(kneeIndex), 'Color', 'r', 'LineWidth', 2);
% Now fit the left portion of the data that follows an exponential to an exponential formula.
coefficients2 = FitExponentialGrowth(T2(1:kneeIndex), L2(1:kneeIndex), 4)
subplot(5, 2, 5);
plot(T3,L3)
grid on;
title('L3 vs. T3', 'FontSize', fontSize);
kneeIndex = FindKneeIndex(L3, T3);
xline(T3(kneeIndex), 'Color', 'r', 'LineWidth', 2);
% Now fit the left portion of the data that follows an exponential to an exponential formula.
coefficients3 = FitExponentialGrowth(T3(1:kneeIndex), L3(1:kneeIndex), 6)
subplot(5, 2, 7);
plot(T4,L4)
grid on;
title('L4 vs. T4', 'FontSize', fontSize);
kneeIndex = FindKneeIndex(L4, T4);
xline(T4(kneeIndex), 'Color', 'r', 'LineWidth', 2);
% Now fit the left portion of the data that follows an exponential to an exponential formula.
coefficients4 = FitExponentialGrowth(T4(1:kneeIndex), L4(1:kneeIndex), 8)
subplot(5, 2, 9);
plot(T5,L5)
grid on;
title('L5 vs. T5', 'FontSize', fontSize);
kneeIndex = FindKneeIndex(L5, T5);
xline(T5(kneeIndex), 'Color', 'r', 'LineWidth', 2);
% Now fit the left portion of the data that follows an exponential to an exponential formula.
coefficients5 = FitExponentialGrowth(T5(1:kneeIndex), L5(1:kneeIndex), 10)
% Maximize figure window.
g = gcf;
g.WindowState = 'maximized';
fprintf('Done running %s.m\n', mfilename);
function kneeIndex = FindKneeIndex(L, T)
slopes = (L(end) - L) ./ (T(end) - T);
kneeIndex = find(slopes < 0.15, 1, 'first');
% figure;
% plot(T, slopes, 'b.-');
% grid on;
end
function coefficients = FitExponentialGrowth(X, Y, plotNumber)
fontSize = 14;
% Convert X and Y into a table, which is the form fitnlm() likes the input data to be in.
tbl = table(X(:), Y(:));
% Define the model as Y = a + exp(-b*x)
% Note how this "x" of modelfun is related to big X and big Y.
% x((:, 1) is actually X and x(:, 2) is actually Y - the first and second columns of the table.
modelfun = @(b,x) b(1) * exp(b(2)*x(:, 1)) + b(3);
beta0 = [11, .5, 0]; % Guess values to start with. Just make your best guess.
% Now the next line is where the actual model computation is done.
mdl = fitnlm(tbl, modelfun, beta0);
% Now the model creation is done and the coefficients have been determined.
% YAY!!!!
% Extract the coefficient values from the the model object.
% The actual coefficients are in the "Estimate" column of the "Coefficients" table that's part of the mode.
coefficients = mdl.Coefficients{:, 'Estimate'}
% Create smoothed/regressed data using the model:
yFitted = coefficients(1) * exp(coefficients(2)*X) + coefficients(3);
% Now we're done and we can plot the smooth model as a red line going through the noisy blue markers.
subplot(5, 2, plotNumber);
plot(X, Y, 'b.', 'MarkerSize', 20);
hold on;
plot(X, yFitted, 'r-', 'LineWidth', 2);
grid on;
title('Exponential Regression with fitnlm()', 'FontSize', fontSize);
xlabel('X', 'FontSize', fontSize);
ylabel('Y', 'FontSize', fontSize);
legendHandle = legend('Noisy Y', 'Fitted Y', 'Location', 'northwest');
legendHandle.FontSize = 12;
formulaString = sprintf('Y = %.3f * exp(%.3f * X) + %.3f', coefficients(1), coefficients(2), coefficients(3))
text(5, 17000, formulaString, 'FontSize', 12, 'FontWeight', 'bold');
end
12 Comments
Mushtaq Al-Jubbori
on 6 Jun 2021
Edited: Mushtaq Al-Jubbori
on 6 Jun 2021
Image Analyst
on 6 Jun 2021
Why??? Both Scott and I suggested a piecewise fit would be more accurate. Sure you could fit a super high order polynomial to it but it would oscillate wildly between the training points - I'm sure you've seen those demos in your classes. So give us some definite rationale for why it's absolutely necessary to have a single formula. You know that a model is not for recreating your training points -- you already know those values -- it's for predicting points that are not in your training set. So why ruin those predictions by having a single formula that won't predict them well.
Mushtaq Al-Jubbori
on 6 Jun 2021
Edited: Mushtaq Al-Jubbori
on 6 Jun 2021
Mushtaq Al-Jubbori
on 6 Jun 2021
Edited: Mushtaq Al-Jubbori
on 6 Jun 2021
Image Analyst
on 6 Jun 2021
Well I think you could put that model into my code for fitnlm() instead of the model I used. Did you do that and try it? How well did it work out.
Mushtaq Al-Jubbori
on 6 Jun 2021
Image Analyst
on 6 Jun 2021
I didn't understand that. So you're going to search for another model? Or you're going to use Alex's model? Please clarify.
Also clarify if you're going to program up your function
"the function A1*tanh(U)^A4, where U=exp((x-A2)/A3)"
Do want to use that equation (instead of another one you're going to search for, or mine, or Alex's model)? If you still want to use your tanh model, let me know if you're going to do it or you need us to do it.
Bottom line, it's not clear if you're going to search for another model, or if you're going to use one mentioned here (either Alex's, mine, Scott's, or your tanh model). What model are you're going to use?
Mushtaq Al-Jubbori
on 7 Jun 2021
Image Analyst
on 7 Jun 2021
Is there a theory that says that your physical process should theoretically follow a coth(x)-1/x curve?
Or is it just something that looks like it fits?
Mushtaq Al-Jubbori
on 7 Jun 2021
Image Analyst
on 8 Jun 2021
Edited: Image Analyst
on 8 Jun 2021
I don't know if coth(x)-1/x is ok or not. It doesn't seem to fit as well as the piecewise exponential growth function I suggested.
x = linspace(-100, 100, 1000);
y = coth(x) - 1 ./ x;
plot(x, y, 'b-', 'lineWidth', 2);
grid on;
but let's say that YOU say it's ok. So please explain the physical meaning of the curve, meaning how it theoretically applies to the physical process you are modeling, which is "Track length of alpha pariclese and etchin time". Presumably someone published a paper where the physics says it should be that equation theoretically.
I don't think the fit looks very good. Defined the model as Y = a * coth(b * x) - c ./ x
Code is attached.
Mushtaq Al-Jubbori
on 8 Jun 2021
Scott MacKenzie
on 6 Jun 2021
Edited: Scott MacKenzie
on 6 Jun 2021
Here's what I put together. This is only for the T5-L5 data. You can repeat this for the rest of the data and build up the plot with the remaining points and lines.
L5=[3 3.4 3.8 4.4 5.2 6 6.8 8 9.2 10.8 12.8 15 16.4 17.4 17.8 17.6 17.8 17.8 17.6 17.8 17.8];%4.85MeV
T5=[2:0.5:7.5 7.7 8:0.5:10.5 11.5 12];
threshold = 0.2; % adjust accordingly
% find elbow (index where curve transitions to flat line)
elbow = find(diff(L5) < threshold, 1);
% plot markers, set axis limits
plot(T5, L5, 'or', 'markerfacecolor', 'r');
ax = gca;
ax.YLim = [0 20];
ax.XLim = [0 15];
% plot flat line
x1 = T5(elbow); x2 = T5(end);
y1 = mean(L5(elbow):L5(end)); y2 = y1;
line([x1 x2], [y1 y2], 'color', 'r', 'linewidth', 1.5);
hold on;
% plot curved line using polynomial fit
x = T5(1:elbow);
ySample = L5(1:elbow);
pf = polyfit(x,ySample,5); % order 5 looks pretty good
y = polyval(pf, x)
y(end) = L5(elbow); % ensure curved line meets flat line
plot(x, y, 'r', 'linewidth', 1.5);
Sulaymon Eshkabilov
on 6 Jun 2021
This nonlinear least squares fit gives quite nice fit models:
...
x = T2;
y = L2;
p = [1, 3]; % Guess: Fit model parameters
q = [-1, -3]; % Guess: Fit model parameters
plot(x,y,'r*')
xlabel 't'
ylabel 'Response'
p = optimvar('p',2);
q = optimvar('q',2);
fun = p(1)*(1+exp(q(1)+q(2)*x) + p(2)*x).^-0.65;
obj = sum((fun - y).^2);
LSQ_Prob = optimproblem("Objective",obj);
x0.p = [1/2, 5/2];
x0.q = [-1/2,-5/2];
show(LSQ_Prob) % Display the problem formulation and parameters
[Sol,fval] = solve(LSQ_Prob,x0) % Show the results
figure
FM_val = evaluate(fun,Sol);
plot(x,y,'r*', x, FM_val,'b-')
legend('Original Data','Fitted Curve', 'location', 'southeast')
xlabel('T')
ylabel('L & Fit Model')
title("Data vs. Fitted Model Response")
SStot = sum((y-mean(y)).^2);
SSres = sum((y-FM_val).^2);
Rsq = 1 - SSres/SStot;
gtext(['R^2 = ' num2str(Rsq)])
gtext(['Fit model: ' num2str(Sol.p(1)) '*(1+exp(' num2str(Sol.q(1)) ' + ' num2str(Sol.q(2)) '*x) + ' num2str(Sol.p(2)) '*x).^{-0.65}'])
fprintf('Found fit model is: %f *(1+exp(%f+%f*x) + %f*x)^-0.65 \n', [Sol.p(1), Sol.q(1), Sol.q(2), Sol.p(2)])
7 Comments
Alex Sha
on 8 Jun 2021
The fitting function of "fun = p(1)*(1+exp(q(1)+q(2)*x) + p(2)*x).^-0.65;" is near same as I post out previously, that is "y = p1/(1+Exp(p2+p3*x)+p4*x)^p5" , unless put p5 as constant of "0.65"
Image Analyst
on 8 Jun 2021
Edited: Image Analyst
on 8 Jun 2021
How did either of you come up with this model? Is it a well known model? Did you just guess? Any theory behind it? Is it a known model for his physical process of "Track length of alpha pariclese and etchin time"?
Alex Sha
on 8 Jun 2021
Hi, the model function was found out automatically by 1stOpt, no any physicall meanning.
Mushtaq Al-Jubbori
on 8 Jun 2021
Image Analyst
on 8 Jun 2021
@Mushtaq Al-Jubbori, I've seen you twice now say that you will explain the theoretical basis for using the coth(x) -1/x function, but you have not done so yet. Why not? Well never mind, I've lost interest. Good luck with your project.
I've given you code to fit that function. Only you can say it it's OK to you and meets your needs.
Mushtaq Al-Jubbori
on 8 Jun 2021
Image Analyst
on 8 Jun 2021
OK, fine. Then use the piecewise function I gave you, or use ad hoc function mysteriously found by that third party software.
but you may need to upgrade your MATLAB, as you said.
Mushtaq Al-Jubbori
on 6 Jun 2021
0 votes
Categories
Find more on Conway's Game of Life in Help Center and File Exchange
on 4 Jun 2021
on 10 Jun 2021
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
