Properly centering interquartile ranges on graphed data

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:

Answer 1

[Turlough Hughes]

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.

[Robert Pearhill]

This is really helpful, thank you so much!