5 views (last 30 days)

Show older comments

I need to fit the attached data as in the image below:

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

Sulaymon Eshkabilov
on 6 Jun 2021

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.

Sulaymon Eshkabilov
on 8 Jun 2021

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

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)])

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.

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

Start Hunting!