How to use fitnlm

Hello I am trying to fit acurve to my data but I don't know where to begin or what to do? Any help will be appretiated!

This is the graph and code of my data. I have attached the excel file of my data.

A=xlsread('Plague_data.xlsx');

t=A(:,12);

c=A(:,13);

figure (1)

Number.of.cases=plot(t,c,'.','color','b')

Answer 1

[Image Analyst]

Try this:

% Fits a Gaussian to the USA daily deaths for CoVid-19, which includes data for only the left portion of the Gaussian, not the full Gaussian.
%--------------------------------------------------------------------------------------------------------------------------------------------------------
% CLEAN UP - INITIALIZATION STEPS
clc;    % Clear the command window.
fprintf('Beginning to run %s.m.\n', mfilename);
close all;  % Close all figures (except those of imtool.)
clearvars;  % Erase all existing variables.
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 20;

%--------------------------------------------------------------------------------------------------------------------------------------------------------
% TRAINING DATA PREPARATION
% Read in data from workbook.
data = readmatrix('Fitnlm_Gaussian_CoVid19.xlsx');

% Get rid of nans
validRows = ~isnan(data(:, 2));
% Get rid of zero counts
validRows = validRows & (data(:, 2) > 0);
data = data(validRows, :)
X = data(:, 1);
Y = data(:, 2);
hFig = figure;
plot(X, Y, 'b.-', 'LineWidth', 2, 'MarkerSize', 30);
grid on;
xlabel('Date', 'FontSize', fontSize);
ylabel('Daily Deaths', 'FontSize', fontSize);
title('CoVid-19 Daily Deaths', 'FontSize', fontSize);
% Make x axis be dates.
datetick('x', 'mmm dd, yyyy', 'keepticks');
ax = gca;
ax.XTickLabelRotation = -45;

% Convert X and Y into a table, which is the form fitnlm() likes the input data to be in.
tbl = table(X(:), Y(:));

%--------------------------------------------------------------------------------------------------------------------------------------------------------
% MODEL DEFINITION
% Define the model as Y = a + b * exp(-(x - c)^2 / d)
% Create an anonymous function for it.
% Note how this "x" of ModelFunction 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.
ModelFunction = @(b, x) b(1) + b(2) * exp(-(x(:, 1) - b(3)).^2/b(4));  

%--------------------------------------------------------------------------------------------------------------------------------------------------------
% MODEL CREATION : ESTIMATION OF PARAMETERS
% Guess starting model values to start with.  Just make your best guess.  
% These are just starting points for the coefficients and will be adjusted during the fit to produce the real coefficients.
beta0 = [0, max(X), mean(X), var(X)]; % A guess at what the coefficients will be.
% Now the next line is where the actual model computation is done.
mdl = fitnlm(tbl, ModelFunction, 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'}

%--------------------------------------------------------------------------------------------------------------------------------------------------------
% MODEL PLOTTING : PLOTTING FITTED/ESTIMATED VALUES OVER THE WHOLE RANGE OF X.
% Let's do a fit, but let's get more points on the fit, beyond just the widely spaced training points,
% so that we'll get a much smoother curve.
xFitted = linspace(min(X), max(X), 1920); % Let's use 1920 points, which will fit across an HDTV screen about one sample per pixel.

% Create smoothed/regressed data using the model:
yFitted = ModelFunction(coefficients, xFitted(:));
% yFitted = coefficients(1) + coefficients(2) * exp(-(xFitted - coefficients(3)).^2 / coefficients(4));
% Now we know that since it's a Gaussian there should not be any negative values.  So clip to 0
yFitted = max(0, yFitted);
% Now we're done and we can plot the smooth model as a red line going through the noisy blue markers.
hold on;
plot(xFitted, yFitted, 'r-', 'LineWidth', 2);
grid on;
title('CoVid-19 Daily Deaths in the USA -- Exponential Regression with fitnlm()', 'FontSize', fontSize);
% Put a line at the peak.
darkGreen = [0, 0.5, 0];
xline(coefficients(3), 'Color', darkGreen, 'LineWidth', 3);
% Put up text for the green line.
str = sprintf(' Peak at %s', datestr(coefficients(3)));
text(coefficients(3), 100, str, 'Color', darkGreen, 'FontSize', 14, 'FontWeight', 'bold');
legendHandle = legend('Actual Y', 'Fitted Y', 'Location', 'northwest');
legendHandle.FontSize = 25;

hFig.WindowState = 'maximized';
fprintf('Done running %s.m.\n', mfilename);

[Peter Derivolkov]

Thank you so much this is perfect!

[Image Analyst]

Attached is a new version where it asks you if you want to allow an offset on the Gaussian or not. Plus it uses a bar chart and has more recent data.

Answer 2

[Keegan Carvalho]

Firstly, you need to know what fitting curve you are looking for; linear, quadratic, etc. I've gone with polyfit which you can later use to evaluate on the variables T and C.

A=xlsread('Plague_data.xlsx'); 
t=A(:,12);
c=A(:,13);
C = isnan(c);
T = isnan(t);                                   % polyfit cannot operate if variables have NaN values
fit1 = polyfit(t(~T),c(~C),1);                  % polyfit can be used and NaN values are omitted

You can start of with the above and later evaluate using polyval. Read about polyfit, fit.

Hope this helps!

[Image Analyst]

If it's a typical pandemic curve, it's probably like the left portion of a Gaussian.

[Peter Derivolkov]

Thank you so much I appreciate it!

[Khadija]

Hi, do you have any idea about fitnlm with function with ODE form?

[Torsten]

You should be able to evaluate the solution structure "mdl" following the documentation.

I think you won't be able to set lower and upper bounds with fitnlm because of the statistical evaluation. But you could square the parameters in the ordinary differential equations.

Another poster in the forum always suggests using "ga" for complicated problems - maybe it's worth trying.

mdl = fitnlm([tdata;tdata],[Hdata;HSdata],@simulatedhs,k0)

function simulated_data = simulatedhs(k,tdata)
  gamma = 1.5;
  phi_S =k(1);
  phi_H =k(2);
  c1 = k(3);
  c2=k(4) ;
  theta1 =k(5) ;
  theta2 = k(6);
  alpha = k(7);
  beta = k(8);
  rho1 = k(9);
  rho2= k(10);
  rho3=k(11);
  [T, Y] = ode15s(@(t,y)modelhs(t,y,k),tdata(1:end/2),[ 5000.0 20.0 2.0  0.0  6.0 2.0  0.0 ]);
  M = Y(:,1);
  S = Y(:,2);
  H1 = Y(:,3);
  H2 = Y(:,4);
  H1S=Y(:,5);
  H2S=Y(:,6);
  H = phi_H * ( H1 + c1 *  H1S ).*M + rho1 * gamma*H1S + alpha* H2; % new cases from mono-HIV
  HS = theta1*phi_S * ( S + c2 *  H1S ).*H1 + theta2*phi_H * ( H1 + c1 *  H1S ).*S+ rho2*alpha*H1S;   %new case of coinfection hiv+syphilis
  simulated_data = [H;HS];
end

[Khadija]

Hi sir,

I hope you are in a good mood, what is the using "ga" for complicated problems?

im trying the last commande and it gives this error:

Error using nlinfit (line 219)

MODELFUN must be a function that returns a vector of fitted values the same size as Y (22-by-1). The

model function you provided returned a result that was 11-by-2.

One common reason for a size mismatch is using matrix operators (*, /, ^) in your function instead of the

corresponding elementwise operators (.*, ./, .^).

Error in NonLinearModel/fitter (line 1168)

nlinfit(X,y,F,b0,opts,wtargs{:},errormodelargs{:});

Error in classreg.regr.FitObject/doFit (line 94)

model = fitter(model);

Error in NonLinearModel.fit (line 1487)

model = doFit(model);

Error in fitnlm (line 99)

model = NonLinearModel.fit(X,varargin{:});

Error in essaie7 (line 73)

mdl = fitnlm([tdata;tdata],[Hdata;HSdata],@simulatedhs,k0);

Any help???

[Torsten]

simulated_data = [H;HS];

This is 22x1, not 11x2. Did you change as written above in "simulatedhs" ?

[Khadija]

Hi Mr. Torsten,

the problem is that the dimensions output of simulatedhs solutions is not the same of Ydata ([H,HS]).

I thought about using the variable 'tmeasure' to have an output of the same dimension. Did I understand the problem correctly?

% Données spécifiques

specific_data = [

2009 2 8;

2010 10 22;

2011 30 45;

2012 111 75;

2013 125 96;

2014 255 192;

2015 379 227;

2016 384 238

2017 360 279;

2018 399 229;

2019 235 128

];

% Utilisez les données spécifiques

tdata = specific_data(:, 1);

Hdata = specific_data(:, 2);

HSdata = specific_data(:, 3);

tforward = 2009:1:2019;

%tmeasure = [1:100:1001]';

% Initial Parameter values

gamma = 1.5;

phi_S =0.0006; % transmission prob

phi_H = 0.000051; % trans proba

c1=3;

c2=1.5;

theta1 =10; % djustment parameters for syph

theta2 =4; %djustment parameters for hi

alpha = 0.6; % progression rate

beta = 0.2; %Complications rate

rho1 =0.4; % adjustment parameters

rho2 =1.5; % adjustment parameters

rho3 = 1.5; % adjustment parameters

k0 = [phi_S phi_H c1 c2 theta1 theta2 alpha beta rho1 rho2 rho3 ];

% fitting the curve and extract statistical information

mdl = fitnlm([tdata;tdata],[Hdata;HSdata],@simulatedhs,k0);

simulated_data = simulatedhs(k,tforward);

H = simulated_data(:,1);

HS = simulated_data(:,2);

figure(3)

plot (tdata.', HSdata,'r*')

hold on

plot(tforward, HS, 'k--','linewidth',2)

legend('H-S coinfection Data', 'Model estimation with fit');

axis([2009 2019 0 1000]);

figure(4)

plot(tforward, H, 'b--','linewidth',2)

hold on

plot (tdata.', Hdata,'ro')

legend('Mono-Hiv Data', 'Model estimation with fit');

axis([2009 2019 0 1000]);

function dy=modelhs(~,y,k)

delta = 0.0001; % Taux de mortalité

delta_S = 0.0005; % Taux de mort de Syphilis.

delta_H = 0.004;

Lambda =4.04 *100;

gamma=0.5;

phi_S =k(1);

phi_H =k(2);

c1 = k(3);

c2=k(4) ;

theta1 =k(5) ;

theta2 = k(6);

alpha = k(7);

beta = k(8);

rho1 = k(9);

rho2= k(10);

rho3=k(11);

dy = zeros(7,1);

%lambda_s=phi_S * ( y(2) + c2 * y(5))

%lambda_H= phi_H * ( y(3) + c1 * y(5) )

dy(1) = Lambda + gamma * y(2) - (phi_S * ( y(2) + c2 * y(5)) + phi_H * (y(3) + c1 * y(5)) + delta ) * y(1) ;%M

dy(2)= phi_S * ( y(2) + c2 * y(5) ) * y(1) - ( gamma + theta2 * phi_H * ( y(3) + c1 * y(5) ) + delta_S ) * y(2) ;%S

dy(3) = phi_H * ( y(3) + c1 * y(5) ) * y(1) + rho1 * gamma * y(5) - (theta1 * phi_S * ( y(2) + c2 * y(5)) + delta + delta_H + alpha) * y(3) ;%H1

dy(4) = alpha * y(3) - (beta + delta + delta_H) * y(4) ;%H2

dy(5) = theta2 * phi_H * ( y(3) + c1 * y(5) ) * y(2) + ( theta1 * phi_S * ( y(2) + c2 * y(5) ) ) * y(3) - ( rho1 * gamma + rho2 * alpha + delta_H + delta_S + delta ) * y(5) ;%H1S

dy(6)= rho2 * alpha * y(5) - ( rho3 * beta + delta_S + delta_H + delta ) * y(6) ;%H2S

dy(7)= beta * y(4) + rho3 * beta * y(6) - ( delta + delta_H ) * y(7) ;%C

end

function simulated_data = simulatedhs(k,tdata)

%tmeasure = [1:100:1001]'; %for evaluating solutions of simulatedhs in tdata

gamma = 1.5;

phi_S =k(1);

phi_H =k(2);

c1 = k(3);

c2=k(4) ;

theta1 =k(5) ;

theta2 = k(6);

alpha = k(7);

beta = k(8);

rho1 = k(9);

rho2= k(10);

rho3=k(11);

[T, Y] = ode15s(@(t,y)modelhs(t,y,k),tdata(1:end/2),[ 5000.0 20.0 2.0 0.0 6.0 2.0 0.0 ],odeset('RelTol',1e-10,'AbsTol',1e-10));

%[T, Y] = ode15s(@(t,y)modelhs(t,y,k),tdata,[ 5000.0 20.0 2.0 0.0 6.0 2.0 0.0 ]);

M = Y(:,1); % M = Y(tmeasure(:),1); something like this!!

S = Y(:,2);

H1 = Y(:,3);

H2 = Y(:,4);

H1S=Y(:,5);

H2S=Y(:,6);

H=phi_H * ( H1 + c1 * H1S ).*M + rho1 * gamma*H1S + alpha* H2; % new cases from mono-HIV

HS=theta1*phi_S * ( S + c2 * H1S ).*H1 + theta2*phi_H * ( H1 + c1 * H1S ).*S+ rho2*alpha*H1S; %new case of coinfection hiv+syphilis

simulated_data = [H,HS];

end

[Torsten]

This is not the function "simulatedhs" I posted above.

Hint: Look at the last line.

Maybe you want to improve your MATLAB skills before continuing with the mathematical challenge. MATLAB offers an online tutorial free-of-costs to learn the basics of the language:

https://matlabacademy.mathworks.com/details/matlab-onramp/gettingstarted

[Khadija]

I aded the semicolon, when running, it gives this error:

Warning: Failure at t=2.009916e+03. Unable to meet integration tolerances without reducing the step size

below the smallest value allowed (3.637979e-12) at time t.

[Torsten]

Yes, that's the same message that I got. "ode15s" is not able to integrate your differential equations beyond t = 2.009916e3 because it cannot cope with the k-vector supplied by "fitnlm".

Note that you cannot set bounds on the parameters in "fitnlm". This means that "fitnlm" might come up with k-vectors that are not physical and lead to difficulties in the integration process.