Hello, I would like to know how to draw a shaded region like this?

I am trying to estimate the parameters for a mathematical model, I managed to do it but there remains a problem of determining the confidence interval and drawing the shaded region which contains the majority of the data points.

[Image Analyst]

If you have any actual questions, then ask them and attach your data and code to read it in with the paperclip icon after you read this:

TUTORIAL: How to ask a question (on Answers) and get a fast answer

In the meantime, see the FAQ: https://matlab.fandom.com/wiki/FAQ#How_do_I_shade_the_region_between_two_curves?

Answer 1

[Star Strider]

Perhaps something like this —

x = linspace(0, 5, 25).';
y = 5*exp(-0.5*x) + randn(size(x));

fcn = @(b,x) b(1).*exp(b(2).*x) + b(3);

B0 = randn(3,1);
mdl = fitnlm(x, y, fcn, B0)

mdl = 
Nonlinear regression model:
    y ~ b1*exp(b2*x) + b3

Estimated Coefficients:
          Estimate      SE        tStat      pValue 
          ________    _______    _______    ________

    b1      5.5856     1.9322     2.8908    0.008483
    b2    -0.30843    0.24258    -1.2715     0.21684
    b3     -1.3477      2.218    -0.6076     0.54968


Number of observations: 25, Error degrees of freedom: 22
Root Mean Squared Error: 0.997
R-Squared: 0.66,  Adjusted R-Squared 0.629
F-statistic vs. constant model: 21.3, p-value = 7.13e-06

[ynew,yci] = predict(mdl, x);

figure
hp1 = plot(x, y, '.', 'DisplayName','Data');
hold on
hp2 = plot(x, ynew, '-r', 'DisplayName','Regression');
hp3 = plot(x, yci, '--r', 'DisplayName','95% CI');
hp4 = patch([x; flip(x)], [yci(:,1); flip(yci(:,2))], 'r', 'FaceAlpha',0.25, 'EdgeColor','none', 'DisplayName','Shaded Confidence Interval');
hold off
grid
xlabel('X')
ylabel('Y')
legend([hp1,hp2,hp3(1),hp4], 'Location','best')

.

[Khadija]

This is almost the result that i am trying to obtain, but i am trying to estimate parameters in differential equations system, the code is

clear all ; clf ; clc ;


% Données spécifiques
specific_data = [
    2013   58041  2030 32755 64;
    2014   52204 2000 36728 68;
    2015   43000  2007 41899 94;
    2016    30400 3124 51385 91;
    2017  24300 3020 55583 97;
    2018   20800 2153  59948 99;
    2019  16613 2000 64288 99;
    2020   13500 2523 65877 100;
    2021   12700 2154 67446 10;
    2022   10150 1473 68279 100
    
    ];

tforward = 2013:1:2022; 
%tmeasure = [1:100:801]'; 
% Utilisez les données spécifiques
tdata = specific_data(:, 1);
Hdata = specific_data(:, 2);
Adata = specific_data(:, 3);
THdata = specific_data(:, 4);
TAdata = specific_data(:, 5);
% Paramètres initiaux

deltaA=0.0003;
delta=0.016+0.0004 ;
deltaC=0.0008 ;
h=1290;
sigma=0.06 ; 
xi= 0.02;  
alpha= 0.04;
epsilon= 0.03;%epsilon=xi*1.5 a revoir la relation
beta= 0.06;% beta=alpha*1.5a revoir la relation
phi= 0.03;%phi=sigma*0.05a revoir la relation
psi= 0.3;
gammaH=0.4 ;
gammaA=0.4 ;
delta=0.0016 ;
deltaA=0.0004 ;
deltaC=0.01 ;
deltaOC= 0.01;
p1= 0.5;
p2=1.5;
p3=2;
k0 = [xi sigma alpha gammaH gammaA deltaOC p1 p2 p3 psi deltaA ];

% Résolution de l'équation différentielle avec les paramètres initiaux
[T, Y] = ode23s(@(t, y) model_CVD(t, y, k0), tforward, [58041  2030 32755 64 1000.0]);
yintH = Y(:,1); 
yintA = Y(:,2);
yintTH = Y(:,3);
yintTA = Y(:,4);


% Affichage des données spécifiques et de la solution de l'équation différentielle
% Display the results

figure(1)
%subplot(1,2,1);
plot(tdata, Hdata, 'r.');
hold on
plot(tdata,  yintH, 'b-');
xlabel('time in days');
ylabel('Number of newly infected cases  ');
legend('Data', 'Model estimation');
axis([2013 2022 0 100000]);
% Minimization routine using Nelder & Mead Simplex algorithm (a derivative-free method)
% Assigns the new values of parameters to k and the error to fval
% Minimization routine using Nelder & Mead Simplex algorithm (a derivative-free method)
% Assigns the new values of parameters to k and the error to fval



[k,fval] = fminsearch(@moderCVD,k0);

%print final values of fitted parametres
disp(k);
xi=k(1);
sigma=k(2);
alpha= k(3);
gammaH=k(4) ;
gammaA=k(5) ;
deltaOC= k(6);
p1= k(7);
p2= k(8);
p3= k(9);
psi= k(10);
deltaA=k(11);


%Draw the data with the final ODE
[T Y] = ode45(@(t,y)(model_CVD(t,y,k)),tforward,[  58041  2030 32755 64 1000.0]);
yintH = Y(:,1);
yintA = Y(:,2);
yintTH = Y(:,3);
yintTA = Y(:,4);
%subplot(1,2,2);

figure(2)
plot(tdata,Hdata,'r.');
hold on
plot(tdata,yintH,'g-');
xlabel('Time in days');
ylabel('Number of newly infected cases ');
legend('Data', 'Model estimation');
axis([2013 2022 0 100000]);
function error_in_data = moderCVD(k)
% Données spécifiques
specific_data = [
    2013   58041  2030 32755 64;
    2014   52204 2000 36728 68;
    2015   43000  2007 41899 94;
    2016    30400 3124 51385 91;
    2017  26300 1320 55583 97;
    2018   20800 2153  59948 99;
    2019  16613 2000 64288 99;
    2020   13500 2523 65877 100;
    2021   12700 2154 67446 10;
    2022   10150 1473 68279 100
    ];

tforward = 2013:1:2022; 
%tmeasure = [1:100:901]'; 
% Utilisez les données spécifiques
tdata = specific_data(:, 1);
Hdata = specific_data(:, 2);
Adata=specific_data(:, 3);
THdata=specific_data(:, 4);
TAdata=specific_data(:, 5);
[T, Y] = ode23s(@(t,y)(model_CVD(t,y,k)),tforward,[ 58041.0  2030.0 32755.0 64.0 1000.0]);

Hq = Y(:,1);
Aq = Y(:,2);
THq = Y(:,3);
TAq = Y(:,4);% assignts the y-coordinates of ...
%the solution at
D=mean(Hdata).^2;
A=(Hq - Hdata).^2;
error_in_data =sum(A)./(9*D);
%% 
end

function dy=model_CVD(~,y,k)
delta=0.016+0.0004 ;
deltaC=0.0008 ;
h=1290;

deltaA=k(11);

xi = k(1);
sigma = k(2);
alpha= k(3);
gammaH= k(4);
gammaA= k(5);
deltaOC= k(6);
p1= k(7);
p2= k(8);
p3= k(9);
psi=k(10);

dy = zeros(5,1);
dy(1)= h-(delta+sigma+xi+gammaH)*y(1);
dy(2)=sigma*y(1) -(delta+alpha+deltaA+gammaA)*y(2);   
dy(3)=gammaH*y(1)+ psi*y(4) -(delta+deltaOC+p1*sigma+p3*xi)*y(3); 
dy(4)=gammaA*y(2)+ p1*sigma*y(3) -(psi+delta+deltaA+p2*alpha)*y(4);
dy(5)=xi*y(1)+alpha*y(2)+p2*alpha*y(4)+p3*xi*y(3) -(delta+deltaC)*y(5);
end

[Star Strider]

In order to calculate the prediction intervals, it is necessary to put ‘model_CVD’ in a form that the Optimization Toolbox function lsqcurvefit function will accept and works with, since I can then use it to calculate the prediction intervals.

Unfortuantely, that does not appear to be possible. One way to make your function work with lsqcurvefit is to put your functions and data in a form that would work with my code in Parameter Estimation for a System of Differential Equations since it uses lsqcurvefit and generally works quite well. You would then use the Statistics and Machine Learning Toolbox function nlpredci to calculate the prediction intervals with tthe information that lsqcurvefit returns.

I worked for a while to get your function to work with lsqcurvefit, however I could not. If you can get your functions and data to work with my code, I can then help you with the prediction intervals.

[Khadija]

I will try to run my function and my data on your code to see the result, hoping that I will succeed, thank you once again.

[Star Strider]

I just now noticed a problem with your system and data. You have 10 data rows (elements), however you are fitting 11 parameters to them. (I added 5 more as initial conditions to your differential equations, however I can deal with those after I estimate them.) The nlpredci function checks that, and refuses to estimate the prediction intervals in a system with more parameters than data elements. We might be able totwork with thtis if you have more data (either before 2013 or after 2022) however not as it currently stands.

TThat aside, I was able to get your system to work with my code (and lsqcurvefit) to get a reasonable fit to your data —

specific_data = [2013   58041  2030 32755 64;
                 2014   52204 2000 36728 68;
                 2015   43000  2007 41899 94;
                 2016    30400 3124 51385 91;
                 2017  24300 3020 55583 97;
                 2018   20800 2153  59948 99;
                 2019  16613 2000 64288 99;
                 2020   13500 2523 65877 100;
                 2021   12700 2154 67446 10;
                 2022   10150 1473 68279 100];

t = specific_data(:,1);
c = specific_data(:,2:end);

% theta0 = [rand(11,1); c(1,:).'];
theta0 = [0.00211024675897065        0.0848635907256016       0.00104703355180015          26.9475035484427      0.174524482995563        0.0157328273883729          5.66807118486959       0.00387034230923033    0.00833057438988443      1.19541412178633e-06        0.0100515768049107          59634.9739810393    981.620191917828          28904.9057150045          1.76409872028741].';

[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c,zeros(size(theta0)));

Solver stopped prematurely.

lsqcurvefit stopped because it exceeded the function evaluation limit,
options.MaxFunctionEvaluations = 1.500000e+03.

fprintf(1,'\tRate Constants:\n')

	Rate Constants:

for k1 = 1:numel(theta)
    fprintf(1, '\t\tk(%2d) = %13.5f\n', k1, theta(k1))
end

		k( 1) =       0.00204
		k( 2) =       0.00636
		k( 3) =       0.00105
		k( 4) =      26.12180
		k( 5) =       0.17441
		k( 6) =       0.01512
		k( 7) =      76.04608
		k( 8) =       0.00375
		k( 9) =       0.00807
		k(10) =       0.00000
		k(11) =       0.01006
		k(12) =   59634.27241
		k(13) =     981.59528
		k(14) =   28910.89007
		k(15) =       1.71120

tv = linspace(min(t), max(t));
Cfit = kinetics(theta, tv);
Cfitt = kinetics(theta,t);

figure
hd = plot(t, c, 'p');
for k1 = 1:size(c,2)
    CV(k1,:) = hd(k1).Color;
    hd(k1).MarkerFaceColor = CV(k1,:);
end
hold on
hlp = plot(tv, Cfit);
for k1 = 1:size(c,2)
    hlp(k1).Color = CV(k1,:);
end
hold off
grid
xlabel('Time')
ylabel('Concentration')
lgdstr = {'H','A','TH','TA'};
legend(hlp, lgdstr, 'Location','best')

figure
tiledlayout(2,2)
for k = 1:size(c,2)
    nexttile
    plot(t, c(:,k), '.')
    hold on
    plot(tv, Cfit(:,k),'-')
    hold off
    grid
    xlabel('Year')
    ylabel('Amplitude')
    title(lgdstr{k})
end

figure
tiledlayout(2,2)
for k2 = 1:size(c,2)
    nexttile
    col = k2+1;
    mdl = fitlm(t,Rsd(:,k2));
    [ynew,yci] = predict(mdl, t);
    hd = plot(t, c(:,k2), 'p');
    % for k1 = 1:size(c,2)
    %     CV(k1,:) = hd(k1).Color;
    %     hd(k1).MarkerFaceColor = CV(k1,:);
    % end
    hold on
    plot(t, ynew+Cfitt(:,k2),'-b', 'LineWidth',1)
    hlp = patch([t; flip(t)], [(Cfitt(:,k2)+yci(:,1)); flip(Cfitt(:,k2)+yci(:,2))], 'b', 'FaceAlpha',0.25, 'EdgeColor','none');
    % for k1 = 1:size(c,2)
    %     hlp(k1).Color = CV(k1,:);
    % end
    hold off
    grid
    xlabel('Time')
    ylabel('Concentration')
end

function C=kinetics(k,t)
    c0=k(end-4:end);
    [T,Cv]=ode45(@(t,y)model_CVD(t,y,k),t,c0);
    
    % function dC=DifEq(t,c)
    %     dcdt=zeros(4,1);
    %     dcdt(1)=-theta(1).*c(1)-theta(2).*c(1);
    %     dcdt(2)= theta(1).*c(1)+theta(4).*c(3)-theta(3).*c(2)-theta(5).*c(2);
    %     dcdt(3)= theta(2).*c(1)+theta(3).*c(2)-theta(4).*c(3)+theta(6).*c(4);
    %     dcdt(4)= theta(5).*c(2)-theta(6).*c(4);
    %     dC=dcdt;
    % end

    function dy=model_CVD(t,y,k)
        delta=0.016+0.0004 ;
        deltaC=0.0008 ;
        h=1290;
        deltaA=k(11);
        xi = k(1);
        sigma = k(2);
        alpha= k(3);
        gammaH= k(4);
        gammaA= k(5);
        deltaOC= k(6);
        p1= k(7);
        p2= k(8);
        p3= k(9);
        psi=k(10);
        dy = zeros(5,1);
        dy(1)= h-(delta+sigma+xi+gammaH)*y(1);
        dy(2)=sigma*y(1) -(delta+alpha+deltaA+gammaA)*y(2);   
        dy(3)=gammaH*y(1)+ psi*y(4) -(delta+deltaOC+p1*sigma+p3*xi)*y(3); 
        dy(4)=gammaA*y(2)+ p1*sigma*y(3) -(psi+delta+deltaA+p2*alpha)*y(4);
        dy(5)=xi*y(1)+alpha*y(2)+p2*alpha*y(4)+p3*xi*y(3) -(delta+deltaC)*y(5);
    end

    C=Cv(:,2:5);
end

% function C=kineticsc(k,t,col)
%     c0=k(end-4:end);
%     [T,Cv]=ode45(@(t,y)model_CVD(t,y,k),t,c0);
% 
%     % function dC=DifEq(t,c)
%     %     dcdt=zeros(4,1);
%     %     dcdt(1)=-theta(1).*c(1)-theta(2).*c(1);
%     %     dcdt(2)= theta(1).*c(1)+theta(4).*c(3)-theta(3).*c(2)-theta(5).*c(2);
%     %     dcdt(3)= theta(2).*c(1)+theta(3).*c(2)-theta(4).*c(3)+theta(6).*c(4);
%     %     dcdt(4)= theta(5).*c(2)-theta(6).*c(4);
%     %     dC=dcdt;
%     % end
% 
%     function dy=model_CVD(t,y,k)
%         delta=0.016+0.0004 ;
%         deltaC=0.0008 ;
%         h=1290;
%         deltaA=k(11);
%         xi = k(1);
%         sigma = k(2);
%         alpha= k(3);
%         gammaH= k(4);
%         gammaA= k(5);
%         deltaOC= k(6);
%         p1= k(7);
%         p2= k(8);
%         p3= k(9);
%         psi=k(10);
%         dy = zeros(5,1);
%         dy(1)= h-(delta+sigma+xi+gammaH)*y(1);
%         dy(2)=sigma*y(1) -(delta+alpha+deltaA+gammaA)*y(2);   
%         dy(3)=gammaH*y(1)+ psi*y(4) -(delta+deltaOC+p1*sigma+p3*xi)*y(3); 
%         dy(4)=gammaA*y(2)+ p1*sigma*y(3) -(psi+delta+deltaA+p2*alpha)*y(4);
%         dy(5)=xi*y(1)+alpha*y(2)+p2*alpha*y(4)+p3*xi*y(3) -(delta+deltaC)*y(5);
%     end
% 
%     C=Cv(:,col);
% end

Tthe equations for ‘H’ and ‘TH’ fit quite well, however your system does not fit ‘TA’ well at all (nor for that matter, ‘A’), so it might be worthwhile checking it to be certain that those are coded correctly. (I have left my other functions in and commented-out so you can experiment with them.)

This is the best I can do. You need more data in order to calculate the prediction intervals. (I have a version of this coded and saved in MATLAB Online for my reference.)

EDIT — (8 Jul 2024 at 13:28)

I cannot get nlpredci to creatte prediction limits for a specific variable. I remember doing it at one point, however I cannot find that Answer, and I apparently did not archive that code. This is probably not possible.

The only option is to just ‘kludge’ it by fitting the residuals to a linear model and then warping the linear model to fit the nonlinear regression.

.

[Khadija]

Hi Mr Sniper, I'm trying to understand your code. One thing I don't quite understand is:

1- the difference between Cfit and Cfitt?

2- the choice of (theta0)?

Thank you for your attention!!

[Star Strider]

Here, ‘Cfit’ is the high-resolution version, and uses ‘tv’ for itts independent variable argument, while ‘Cfitt’ uses the original time argument, ‘t’.

tv = linspace(min(t), max(t));

Cfit = kinetics(theta, tv);

Cfitt = kinetics(theta,t);

The reason for ‘Cfitt’ is that it matches the original dimensions used in the prediction interval calculations. They use the ‘Rsd’ residuals matrix, since that is the only way I could make this work, using fittlm with ‘t’ and the ‘Rsd’ columns, and predict to generate the prediction intervals. (Prediction intervals are not the same as confidence intervals, however MATLAB only calculates predicttion intervals. I intend to finish writing a File Excnage conttributton that will take tthe prediction intervals and calculate confidence intervals from them, however I’m not there yet.) As I mentioned, calculating them this way is a ‘kludge’, however it is the only way I could calculate the prediction intervals, since it is not possible to generate prediction intervals from a matrix lsqcurvefit output. There is no MATLAB function for that. Plotting the patch objects (the filled prediction intervals) was straightforward, once I figured out the rest of it.

The original ‘theta0’ was random. When I got parameters that worked (and gave good results), I used them simply for efficiency. It takes a few seconds tot run the code with the specific parameters, and much longer with a random vector thtat does not always converge on the same parameter estimates. I wanted to get the code running without waiting an unreasonably long time for it to converge, and I wanted it to converge on the same working parameters in each run.

.

[Khadija]

Hi,

I only have 11 parameters to estimate. While theta0 contains 15 elements

[Star Strider]

First, the problem is that you are estimating 11 parameters witth 10 data sets, The resulting parameters cannot be uniquely estimated as the result. (Ideally, tthere should be at least 1 more data set than parametters estimated.)

My code treats the initial conditions for the differential equations as separate parameters:

c0=k(end-4:end);

(specifically the last 5) and estimates them as such. (They are not kinetic parameters. Treattting the initial conditions as parameters makes it easier to fit the differenttial equations.)

Second, while you have 5 differential equations, you can only fit 4 of them since you only have 4 dependent variables, those being:

specific_data(:,2:5);

The time vector is the first column, and cannot be fit since it is the independent variable. After experimenting a bit, I chose the 4 that provided the best fit to the data.

Your code has problems that I have done my best to solve. You need more data or fewer parameters if you want to estimate the parameters appropriately and reliably and have them actually mean something.

I can go no further on this. I have done everything you requested.

.

[Khadija]

Thank you very much for your patience!!

[Star Strider]

As always, my pleasure!