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

Question

Khadija 2024-7-5，18:07

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2134866-hello-i-would-like-to-know-how-to-draw-a-shaded-region-like-this

评论： Star Strider 2024-7-14，1:02

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.

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Image Analyst 2024-7-5，18:33

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 2024-7-5，18:50

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2134866-hello-i-would-like-to-know-how-to-draw-a-shaded-region-like-this#answer_1481816

在 MATLAB Online 中打开

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

.

10 个评论
显示 8更早的评论隐藏 8更早的评论

Khadija 2024-7-5，19:28

在 MATLAB Online 中打开

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 2024-7-7，18:56

编辑：Star Strider 2024-7-8，13:27

在 MATLAB Online 中打开

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.

.

Star Strider 2024-7-13，1:45

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.

.

Star Strider 2024-7-13，23:59

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 2024-7-14，0:41