Properly centering interquartile ranges on graphed data

18 views (last 30 days)
Hi All,
I'm back again with another question. I am attempting to graph the median set of values from 1,000 simulations of an ODE with accompanying error bars representing the interquartile range. I have been able to succesfully graph the median and generate the interquartile ranges, but using the "errorbar" function as shown below takes my interquartile ranges and centers them symmetrically about the median such that they end up reaching into the negative range. Done correctly, I imagine the error bars should look something like the image I've attached to this post. Thanks for all of the help, this site is truly the bomb!
clc
clearvars
%Total number of replicates
nr=1000;
%Time step and length
tspan=[(1:1:270)];
samples = length(tspan);
result=zeros(samples,nr); % simple numerical 2D array ; preallocation
R0_result=zeros(nr,1);
%Figure hold to plot all 1000 replicates on a single graph
figure; hold on
n=16956; %
y20=0;
tr=0.0000001;
for k=1:nr
%Parameters
y30=(50-20)*rand+20;
p=(407-341)*rand+341;
d=(1/42-1/50)*rand+1/50;
yc=(4.923*10^-3-8.8*10^-5)*rand+8.8*10^-5;
yi=(4.923*10^-3-8.8*10^-5)*rand+8.8*10^-5;
a=(1/30-1/120)*rand+1/120;
r=(0.30-.10)*rand+0.1;
i=(1/4-1/15)*rand+1/15;
bc=(1.5*10^-3-1.5*10^-5)*rand+1.5*10^-5;
bi=(1.5*10^-3-1.5*10^-5)*rand+1.5*10^-5;
c=(50-5)*rand+5;
%Function definition
[t,y] = ode23s(@(t,y) [p*(1-yc-yi) + a*y(2) + tr*y(3) - (c*bc*y(2)*y(1))/n - (c*bi*y(3)*y(1))/n - d*y(1); p*yc + (c*bc*y(2)*y(1))/n + (c*bi*y(3)*y(1))/n - a*y(2) - r*i*y(2) - d*y(2); p*yi + r*i*y(2) + - d*y(3)-tr*y(3);], tspan, [n y20 y30]);
result(:,k) = y(:,3); % simple num array
end
%Translate double array into timeseries for use in the iqr function
Tresult=timeseries(result,length(tspan));
y3_median = median(result,2);
y3_iqr=iqr(Tresult);
errorbar(t,y3_median,y3_iqr, 'bp')
plot(t, y3_median,'r','Linewidth',4);
ylim([0 200])
xlim([0 300])
hold off
Example Graph:

Accepted Answer

Turlough Hughes
Turlough Hughes on 8 Dec 2021
You can use the quantile function
Your simulation:
clc
clearvars
%Total number of replicates
nr=1000;
%Time step and length
tspan=[(1:1:270)];
samples = length(tspan);
result=zeros(samples,nr); % simple numerical 2D array ; preallocation
R0_result=zeros(nr,1);
n=16956; %
y20=0;
tr=0.0000001;
for k=1:nr
%Parameters
y30=(50-20)*rand+20;
p=(407-341)*rand+341;
d=(1/42-1/50)*rand+1/50;
yc=(4.923*10^-3-8.8*10^-5)*rand+8.8*10^-5;
yi=(4.923*10^-3-8.8*10^-5)*rand+8.8*10^-5;
a=(1/30-1/120)*rand+1/120;
r=(0.30-.10)*rand+0.1;
i=(1/4-1/15)*rand+1/15;
bc=(1.5*10^-3-1.5*10^-5)*rand+1.5*10^-5;
bi=(1.5*10^-3-1.5*10^-5)*rand+1.5*10^-5;
c=(50-5)*rand+5;
%Function definition
[t,y] = ode23s(@(t,y) [p*(1-yc-yi) + a*y(2) + tr*y(3) - (c*bc*y(2)*y(1))/n - (c*bi*y(3)*y(1))/n - d*y(1); p*yc + (c*bc*y(2)*y(1))/n + (c*bi*y(3)*y(1))/n - a*y(2) - r*i*y(2) - d*y(2); p*yi + r*i*y(2) + - d*y(3)-tr*y(3);], tspan, [n y20 y30]);
result(:,k) = y(:,3); % simple num array
end
%Translate double array into timeseries for use in the iqr function
Tresult=timeseries(result,length(tspan));
y3_median = median(result,2);
Plotting upper and lower quartiles with errorbar:
Yq = quantile(result,[.25 .75],2);
figure(), errorbar(t, y3_median, y3_median-Yq(:,1) ,Yq(:,2) - y3_median)
hold on, plot(t, y3_median,'r','Linewidth',4);
If you want something more similar to the graph you attached I suggest either using less samples from your simulation or alternatively using lines for your upper and lower bounds and a filled region in between. Here's an example of the former option:
sample_gap = 15;
t_s = t(1:sample_gap:end);
y3_median_s = y3_median(1:sample_gap:end);
Yq_s = Yq(1:sample_gap:end,:);
hf = figure();
errorbar(t_s, y3_median_s,y3_median_s-Yq_s(:,1) ,Yq_s(:,2)-y3_median_s,...
'ok', 'MarkerFaceColor', [0.8500 0.3250 0.0980], 'LineWidth', 2)
set(gca, 'color', [0.8275 0.8980 0.9098])
hf.Position(4) = hf.Position(4)-150;
box off
Note, the quality of this graph in the browser is reduced... the quality in Matlab will be much better.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!

Translated by