# Properly centering interquartile ranges on graphed data

22 views (last 30 days)
Robert Pearhill on 8 Dec 2021
Commented: Robert Pearhill on 8 Dec 2021
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: Turlough Hughes on 8 Dec 2021
You can use the quantile function
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.
Robert Pearhill on 8 Dec 2021
This is really helpful, thank you so much!