Code covered by the BSD License

# Violin Plots for plotting multiple distributions (distributionPlot.m)

by

### Jonas (view profile)

13 Apr 2009 (Updated )

Function for plotting multiple histograms side-by-side in 2D - better than boxplot.

### Editor's Notes:

This file was selected as MATLAB Central Pick of the Week

distributionPlot(varargin)
```function handles = distributionPlot(varargin)
%DISTRIBUTIONPLOT creates violin plots for convenient visualization of multiple distributions
%
% SYNOPSIS: handles = distributionPlot(data,propertyName,propertyValue,...)
%           handles = distributionPlot(ah,...)
%
% INPUT data : m-by-nData array of values, or vector of grouped data (use
%           the 'groups' property to specify the grouping variable), or
%           cell array of length nData.
%           The cell array can either contain vectors with values, or
%           m-by-2 arrays with [bins,counts] if you want to determine the
%           histograms by yourself (m can be different between cell
%           elements). Note that arrays inside cells with any
%           other shape than m-by-2 are reshaped to vector an a warning is
%           thrown (DISTRIBUTIONPLOT:AUTORESHAPE).
%
%       DISTRIBUTIONPLOT accepts the following propertyName/propertyValue
%           pairs (all are optional):
%
%       distWidth :  width of distributions; ideally between 0 and 1.
%           1 means that adjacent distributions might touch. Default: 0.9
%       variableWidth : If true, the width of the distribution changes,
%           reflecting the shape of the histogram of the data. If false,
%           the distribution is only encoded by color levels. Default: true
%       color : uniform coloring of histograms. Supply either a color
%           string ('r'), or a truecolor vector ([1 0 0]). Use a
%           cell array of length nData to specify one color per
%           distribution. Default: 'k'
%           If variableWidth is set to false, a colormap is generated that
%           goes from white to the chose color (or from black, if
%           invert==true).
%           If both 'color', and 'colormap' are specified, 'colormap' takes
%           precedence.
%       colormap : colormap used to describe the distribution (first row
%           corresponds to bins with least data, last row corresponds to
%           bins with most data (invert the grayscale colormap to have
%           black indicate the most data).
%           Supply a cell array of length nData to color distributions
%           individually. Note that using multiple colormaps means that
%           the colorbar doesn't contain much useful information.
%           Default: []
%           Colormap will index into the figure colormap, which will be
%           modified by distributionPlot. This is done to allow editing the
%           distributions in e.g. Adobe Illustrator.
%           If both 'color', and 'colormap' are specified, 'colormap' takes
%           precedence.
%       globalNorm : normalization for bin width (x-direction)
%           0 : every histogram is normalized individually so that the
%               maximum bin width is equal to distWidth. This is best
%               suited to comparing distribution shapes. Default.
%           1 : histograms are normalized such that equal bin width
%               reports equal numbers of counts per bin.
%           2 : histograms are normalized so that the relative areas
%               covered by the histograms reflect the relative total number
%               of data points.
%           3 : histograms areas are normalized so that relative densities
%               are the same across histograms. Thus, if
%               data = {rand(100,1),rand(500,1)},
%               then
%               distributionPlot(data,'globalNorm',2,'histOpt',0,'divFactor',10)
%               shows the left histogram 5x as wide as the right, while
%               distributionPlot(data,'globalNorm',3,'histOpt',0,'divFactor',10)
%               displays both histograms equally wide, since each bin
%               contains ~10% of the data.
%           Options 1 and 2 produce similar results if the bins are spaced
%           equally for the distributions. Options 0 and 3 produce similar
%           results if the data are drawn from the same distributions.
%           Note that colormaps currently always report the number of data
%           points per bin; 'globalNorm' only applies to the distribution
%           shape.
%
%       groups : grouping variable for grouped data. Grouping will be
%                   resolved by calling grp2idx, and unless xNames have
%                   been supplied, group names determine the x-labels.
%                   If the grouping variable is numeric, group labels also
%                   determine x-values, unless the parameter xValues has
%                   been specified.
%       histOpt : histogram type to plot
%                   0 : use hist command (no smoothing, fixed number of
%                       bins)
%                   1 : smoothened histogram using ksdensity with
%                       Normal kernel. Default.
%                   1.1: smoothened histogram using ksdensity where the
%                       kernel is robustly estimated via histogram.m.
%                       Normal kernel.
%                   2 : histogram command (no smoothing, automatic
%                       determination of thickness (y-direction) of bins)
%       divFactor : Parameter dependent on histOpt. If...
%                   histOpt == 0: divFactor = # of bins. Default: 25.
%                       Alternatively, pass a vector which will be
%                       interpreted as bin centers.
%                   histOpt == 1: divFactor decides by how much the default
%                       kernel-width is multiplied in order to avoid an
%                       overly smooth histogram. Default: 1/2
%                   histOpt == 2: divFactor decides by how much the
%                       automatic bin width is multiplied in order to have
%                       more (<1) or less (>1) detail. Default: 1
%       addSpread : if 1, data points are plotted with plotSpread.
%                   distWidth is ideally set to 0.95
%                   This option is not available if the data is supplied as
%                   histograms.
%                   Exchange using the link in the remarks
%       showMM : if 1, mean and median are shown as red crosses and
%                green squares, respectively. This is the default
%                2: only mean
%                3: only median
%                4: mean +/- standard error of the mean (no median)
%                5: mean +/- standard deviation (no median)
%                6: draw lines at the 25,50,75 percentiles (no mean)
%                0: plot neither mean nor median
%       xValues: x-coordinate where the data should be plotted.
%                If xValues are given, "distWidth" is scaled by the median
%                difference between adjacent (sorted) x-values. Note that
%                this may lead to overlapping distributions. Default:
%                1:nData
%       xNames : cell array of length nData containing x-tick names
%               (instead of the default '1,2,3')
%       xMode  : if 'auto', x-ticks are spaced automatically. If 'manual',
%                there is a tick for each distribution. If xNames is
%                provided as input, xMode is forced to 'manual'. Default:
%                'manual'.
%          NOTE: SPECIFYING XNAMES OR XVALUES OR XMODE WILL ERASE PREVIOUS
%                LABELS IF PLOTTING INTO EXISTING AXES
%       yLabel : string with label for y-axis. Default : ''
%                If empty and data is histograms, ylabel is set to 'counts'
%       invert : if 1, axes color is changed to black, and colormap is
%                   inverted.
%       histOri: Orientation of histogram. Either 'center', 'left', or
%                'right'. With 'left' or 'right', the left or right half of
%                the standard violin plot is shown. Has no effect if
%                variableWidth is false. Default: center
%       xyOri  : orientation of axes. Either 'normal' (=default), or
%                'flipped'. If 'flipped', the x-and y-axes are switched, so
%                that violin plots are horizontal. Consequently,
%                axes-specific properties, such as 'yLabel' are applied to
%                the other axis.
%       widthDiv : 1-by-2 array with [numberOfDivisions,currentDivision]
%                widthDiv allows cutting the stripe dedicated to a single
%                distribution into multible bands, which can be filled with
%                sequential calls to distributionPlot. This is one way
%                to compare two (or more) sequences of distributions. See
%                example below.
%       ah : axes handle to plot the distributions. Default: gca
%
% OUTPUT handles : 1-by-4 cell array with patch-handles for the
%                  distributions, plot handles for mean/median, the
%                  axes handle, and the plotSpread-points handle
%
%
% EXAMPLES
%         %--Distributions contain more information than boxplot can capture
%         r = rand(1000,1);
%         rn = randn(1000,1)*0.38+0.5;
%         rn2 = [randn(500,1)*0.1+0.27;randn(500,1)*0.1+0.73];
%         rn2=min(rn2,1);rn2=max(rn2,0);
%         figure
%         ah(1)=subplot(3,4,1:2);
%         boxplot([r,rn,rn2])
%         ah(2)=subplot(3,4,3:4);
%         distributionPlot([r,rn,rn2],'histOpt',2); % histOpt=2 works better for uniform distributions than the default
%         set(ah,'ylim',[-1 2])
%
%         %--- additional options
%
%         data = [randn(100,1);randn(50,1)+4;randn(25,1)+8];
%         subplot(3,4,5)
%
%         %--- defaults
%         distributionPlot(data);
%         subplot(3,4,6)
%
%         %--- show density via custom colormap only, show mean/std,
%         distributionPlot(data,'colormap',copper,'showMM',5,'variableWidth',false)
%         subplot(3,4,7:8)
%
%         %--- auto-binwidth depends on # of datapoints; for small n, plotting the data is useful
%         % note that this option requires the additional installation
%         % of plotSpread from the File Exchange (link below)
%
%         %--- show quantiles
%         subplot(3,4,9),distributionPlot(randn(100,1),'showMM',6)
%
%         %--- horizontal orientation
%         subplot(3,4,10:11),
%         distributionPlot({chi2rnd(3,1000,1),chi2rnd(5,1000,1)},'xyOri','flipped','histOri','right','showMM',0),
%         xlim([-3 13])
%
%         %--- compare distributions side-by-side (see also example below)
%         % plotting into specified axes will throw a warning that you can
%         % turn off using " warning off DISTRIBUTIONPLOT:ERASINGLABELS "
%         ah = subplot(3,4,12);
%         subplot(3,4,12),distributionPlot(chi2rnd(3,1000,1),'histOri','right','color','r','widthDiv',[2 2],'showMM',0)
%         subplot(3,4,12),distributionPlot(chi2rnd(5,1000,1),'histOri','left','color','b','widthDiv',[2 1],'showMM',0)
%
%         %--Use globalNorm to generate meaningful colorbar
%         data = {randn(100,1),randn(500,1)};
%         figure
%         distributionPlot(data,'globalNorm',true,'colormap',1-gray(64),'histOpt',0,'divFactor',[-5:0.5:5])
%         colorbar
%
%         %--Use widthDiv to compare two series of distributions
%         data1 = randn(500,5);
%         data2 = bsxfun(@plus,randn(500,5),0:0.1:0.4);
%         figure
%         distributionPlot(data1,'widthDiv',[2 1],'histOri','left','color','b','showMM',4)
%         distributionPlot(gca,data2,'widthDiv',[2 2],'histOri','right','color','k','showMM',4)
%
%         %--Christmas trees!
%           x=meshgrid(1:10,1:10);
%           xx = tril(x);
%           xx = xx(xx>0);
%           figure
%           set(hh{4}{1},'color','r','marker','o')
% END
%
% REMARKS To show distributions as clouds of points (~beeswarm plot),
%         additional function plotSpread.m from the File Exchange
%
%         I used to run ksdensity with the Epanechnikov kernel. However,
%         for integer data, the shape of the kernel can produce peaks
%         between the integers, which is not ideal (use histOpt=2 for
%         integer valued data).
%
%         A previous iteration of distributionPlot used the input
%         specifications below. They still work to ensure backward
%         compatibility, but are no longer supported or updated.
%           where distWidth of 1 means that the maxima
%           of  two adjacent distributions might touch. Negative numbers
%           indicate that the distributions should have constant width, i.e
%           the density is only expressed through greylevels.
%           Values between 1 and 2 are like values between 0 and 1, except
%           that densities are not expressed via graylevels. Default: 1.9
%
%
%

% created with MATLAB ver.: 7.6.0.324 (R2008a) on Windows_NT
%
% created by: Jonas Dorn; jonas.dorn@gmail.com
% DATE: 08-Jul-2008
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%====================================
%% TEST INPUT
%====================================

% set defaults
def.xNames = [];
def.showMM = 1;
def.distWidth = 0.9;
def.histOpt = 1;
def.divFactor = [25,2,1];
def.invert = false;
def.colormap = [];
def.color = 'k';
def.globalNorm = false;
def.variableWidth = true;
def.groups = [];
def.yLabel = '';
def.xValues = '';
def.xMode = 'manual';
def.histOri = 'center';
def.xyOri = 'normal';
def.widthDiv = [1 1];
isHistogram = false; %# this parameter is not set by input

if nargin == 0 || isempty(varargin{1})
error('not enough input arguments')
end

% check for axes handle
if ~iscell(varargin{1}) && isscalar(varargin{1}) == 1 && ...
ishandle(varargin{1}) && strcmp(get(varargin{1},'Type'),'axes')
ah = varargin{1};
data = varargin{2};
varargin(1:2) = [];
newAx = false;

else
ah = gca;
data = varargin{1};
varargin(1) = [];
newAx = true;
end

% check for current axes limits. Set NaN if the axes have no children
% yet - we need that in case we're building a complicated set of
% distributions
if ~isempty(get(ah,'children'))
xAxLim = xlim;
yAxLim = ylim;
else
[xAxLim,yAxLim] = deal([NaN NaN]);
end

fh = get(ah,'Parent');

% check data. If not cell, convert
if ~iscell(data)
[nPoints,nData] = size(data);
data = mat2cell(data,nPoints,ones(nData,1));
else
% get nData
data = data(:);
nData = length(data);
% make sure all are vectors
badCol = ~cellfun(@isvector,data) & ~cellfun(@isempty,data);
if all(nCols==2)
% bins,counts
isHistogram = true;
else
warning('DISTRIBUTIONPLOT:AUTORESHAPE',...
'Elements %s of the cell array are not vectors. They will be reshaped automatically',...
end
end
end

parserObj = inputParser;
parserObj.FunctionName = 'distributionPlot';
stdWidth = 1; % scaling parameter for variableWidth with uneven x-values
% check whether we're dealing with pN/pV or straight arguments
if ~isempty(varargin) && ~ischar(varargin{1}) && ~isstruct(varargin{1})
% use old format
def.distWidth = 1.9;

parserObj.parse(varargin{:});
opt = parserObj.Results;
% fill in defaults that are not supported in the old version of the
% code
opt.colormap = [];
opt.variableWidth = true;
opt.histOri = 'center';
opt.xValues = [];
opt.xMode = 'auto';
opt.xyOri = 'normal';
opt.widthDiv = [1 1];

% overwrite empties with defaults - inputParser considers empty to be a
% valid input.
fnList = fieldnames(opt);
for fn = fnList'
if isempty(opt.(fn{1}))
opt.(fn{1}) = def.(fn{1});
end
end

% fix a few parameters
if opt.distWidth > 1
opt.distWidth = opt.distWidth - 1;
else
opt.colormap = 1-gray(128);
end
if opt.distWidth < 0
opt.variableWidth = false;
opt.distWidth = abs(opt.distWidth);
end

if ~isempty(opt.xNames)
opt.xMode = 'manual';
end

else
defNames = fieldnames(def);
for dn = defNames(:)'
end

parserObj.parse(varargin{:});
opt = parserObj.Results;

% if groups: deal with data
if ~isempty(opt.groups)
[idx,labels,vals] = grp2idx(opt.groups);
% convert data to cell array
data = accumarray(idx,data{1},[],@(x){x});
nData = length(data);
% if not otherwise provided, use group labels for xnames
if isempty(opt.xNames)
opt.xNames = labels;
if ~iscell(opt.xNames)
opt.xNames = num2cell(opt.xNames);
end
end
if isnumeric(vals) && isempty(opt.xValues)
opt.xValues = vals;
end

end

if ~ischar(opt.xyOri) || ~any(ismember(opt.xyOri,{'normal','flipped'}))
error('option xyOri must be either ''normal'' or ''flipped'' (is ''%s'')',opt.xyOri);
end

end
% common checks

% default x-values: 1:n
if isempty(opt.xValues)
opt.xValues = 1:nData;
elseif length(opt.xValues) ~= nData
error('please supply as many x-data values as there are data entries')
elseif length(opt.xValues) > 1 % only check for scale if more than 1 value
% scale width
stdWidth = median(diff(sort(opt.xValues)));
opt.distWidth = opt.distWidth * stdWidth;
end

if ~isscalar(opt.divFactor) && length(opt.divFactor) == 3 && all(opt.divFactor==def.divFactor)
opt.divFactor = opt.divFactor(floor(opt.histOpt)+1);
end
if isHistogram
opt.histOpt = 99;
if isempty(opt.yLabel)
opt.yLabel = 'counts';
end
end

% check colors/colormaps: do we need to expand colormap?
if ~iscell(opt.colormap)
opt.colormap = {opt.colormap};
end
if ~iscell(opt.color)
opt.color = {opt.color};
end
for iColor = 1:length(opt.color)
if ischar(opt.color{iColor})
opt.color{iColor} = colorCode2rgb(opt.color{iColor});
end
end

% expand - if only single colormap specified, we expand only once
if ~opt.variableWidth
missingColormaps = find(cellfun(@isempty,opt.colormap));
for iMissing = missingColormaps(:)'

endColor = opt.color{max(iMissing,length(opt.color))};
% normally, we go from white to color
cmap = zeros(128,3);
for rgb = 1:3
cmap(:,rgb) = linspace(1,endColor(rgb),128);
end
opt.colormap{iMissing} = cmap;

end
end

% if we have colormaps, we need to create a master which we add to the
% figure. Invert if necessary, and expand the cell array to nData
colormapLength = cellfun(@(x)size(x,1),opt.colormap);
if any(colormapLength>0)

colormap = cat(1,opt.colormap{:});
if opt.invert
colormap = 1-colormap;
end
set(fh,'Colormap',colormap)
if length(opt.colormap) == 1
opt.colormap = repmat(opt.colormap,nData,1);
colormapLength = repmat(colormapLength,nData,1);
colormapOffset = zeros(nData,1);
singleMap = true;
else
colormapOffset = [0;cumsum(colormapLength(1:end-1))];
singleMap = false;
end

else

colormapLength = zeros(nData,1);
if length(opt.color) == 1
opt.color = repmat(opt.color,nData,1);
end
if opt.invert
opt.color = cellfun(@(x)1-x,opt.color,'uniformOutput',false);
end
end

% set hold on
holdState = get(ah,'NextPlot');

% if new axes: invert
if newAx && opt.invert
set(ah,'Color','k')
end

%===================================

%===================================
%% PLOT DISTRIBUTIONS
%===================================

% assign output
hh = NaN(nData,1);
[m,md,sem,sd] = deal(nan(nData,1));
if opt.showMM == 6
md = nan(nData,3,3); % md/q1/q3, third dim is y/xmin/xmax
end

% get base x-array
% widthDiv is a 1-by-2 array with
% #ofDivs, whichDiv
% The full width (distWidth) is split into
% #ofDivs; whichDiv says which "stripe" is active
xWidth = opt.distWidth/opt.widthDiv(1);
xMin = -opt.distWidth/2;
xLow = xMin + xWidth * (opt.widthDiv(2)-1);
xBase = [-xWidth;xWidth;xWidth;-xWidth]/2;
xOffset = xLow + xWidth/2;

% b/c of global norm: loop twice
plotData = cell(nData,2);

% loop through data. Prepare patch input, then draw patch into gca
for iData = 1:nData
currentData = data{iData};
% only plot if there is some finite data
if ~isempty(currentData(:)) && any(isfinite(currentData(:)))

switch floor(opt.histOpt)
case 0
% use hist
[xHist,yHist] = hist(currentData,opt.divFactor);

case 1
% use ksdensity

if opt.histOpt == 1.1
% use histogram to estimate kernel
[dummy,x] = histogram(currentData); %#ok<ASGLU>
if length(x) == 1
% only one value. Make fixed distribution
dx = 0.1;
yHist = x;
xHist = sum(isfinite(currentData));
else
dx = x(2) - x(1);

% make sure we sample frequently enough
x = min(x)-dx:dx/3:max(x)+dx;
[xHist,yHist] = ksdensity(currentData,x,'kernel','normal','width',dx/(1.5*opt.divFactor));
end
else

% x,y are switched relative to normal histogram
[xHist,yHist,u] = ksdensity(currentData,'kernel','normal');
% take smaller kernel to avoid over-smoothing
if opt.divFactor ~= 1
[xHist,yHist] = ksdensity(currentData,'kernel','normal','width',u/opt.divFactor);
end
end

% modify histogram such that the sum of bins (not the
% integral under the curve!) equals the total number of
% observations, in order to be comparable to hist
xHist = xHist/sum(xHist)*sum(isfinite(currentData));

case 2
% use histogram - bar heights are counts as in hist
[xHist,yHist] = histogram(currentData,opt.divFactor,0);
case 99
% bins,counts already supplied
xHist = currentData(:,2)';
yHist = currentData(:,1)';
end
plotData{iData,1} = xHist;
plotData{iData,2} = yHist;
end
end

goodData = find(~cellfun(@isempty,plotData(:,1)));
% get norm
switch opt.globalNorm
case 3
% #3 normalizes relative densities
xNorm(goodData) = cellfun(@(x)min(diff(x)),plotData(goodData,2));
xNorm(goodData) = xNorm(goodData) .* cellfun(@sum,plotData(goodData,1))';
maxNorm(goodData) = cellfun(@max,plotData(goodData,1));
xNorm(goodData) = xNorm(goodData)*max(maxNorm(goodData)./xNorm(goodData));

case 2
% #2 should normalize so that the integral of the
% different histograms (i.e. area covered) scale with the
% respective sum of counts across all bins. Requires evenly spaced
% histograms at the moment
xNorm(goodData) = cellfun(@(x)min(diff(x)),plotData(goodData,2));
maxNorm(goodData) = cellfun(@max,plotData(goodData,1));
xNorm(goodData) = xNorm(goodData)*max(maxNorm(goodData)./xNorm(goodData));
case 1
xNorm(goodData) = max(cat(2,plotData{:,1}));
case 0
xNorm(goodData) = cellfun(@max,plotData(goodData,1));
end

for iData = goodData'

% find current data again
currentData = data{iData};

xHist = plotData{iData,1};
yHist = plotData{iData,2};

% find y-step
dy = min(diff(yHist));
if isempty(dy)
dy = 0.1;
end

% create x,y arrays
nPoints = length(xHist);
xArray = repmat(xBase,1,nPoints);
yArray = repmat([-0.5;-0.5;0.5;0.5],1,nPoints);

% x is iData +/- almost 0.5, multiplied with the height of the
% histogram
if opt.variableWidth

tmp = xArray.*repmat(xHist,4,1)./xNorm(iData);

switch opt.histOri
case 'center'
% we can simply use xArray
xArray = tmp;
case 'right'
% shift everything to the left
delta = tmp(1,:) - xArray(1,:);
xArray = bsxfun(@minus,tmp,delta);
case 'left'
% shift everything to the right
delta = tmp(1,:) - xArray(1,:);
xArray = bsxfun(@plus,tmp,delta);
end

xArray = xArray + opt.xValues(iData);

else
xArray = xArray + iData;
end

% add offset (in case we have multiple widthDiv)
xArray = xArray + xOffset;

% yData is simply the bin locations
yArray = repmat(yHist,4,1) + dy*yArray;

vertices = [xArray(:),yArray(:)];
faces = reshape(1:numel(yArray),4,[])';

if colormapLength(iData) == 0
colorOpt = {'FaceColor',opt.color{iData}};
else
% calculate index into colormap
if singleMap
% use scaled mapping so that colorbar is meaningful
if opt.globalNorm > 0
colorOpt = {'FaceVertexCData',xHist','CDataMapping','scaled','FaceColor','flat'};
else
colorOpt = {'FaceVertexCData',xHist'/xNorm(iData),'CDataMapping','scaled','FaceColor','flat'};
end

else
idx = round((xHist/xNorm(iData))*(colormapLength(iData)-1))+1;
colorOpt = {'FaceVertexCData',idx'+colormapOffset(iData),'CDataMapping','direct','FaceColor','flat'};
end
end

switch opt.xyOri
case 'normal'
hh(iData)= patch('Vertices',vertices,'Faces',faces,'Parent',ah,colorOpt{:},'EdgeColor','none');
case 'flipped'
hh(iData)= patch('Vertices',vertices(:,[2,1]),'Faces',faces,'Parent',ah,colorOpt{:},'EdgeColor','none');
end

if opt.showMM > 0
if isHistogram
[m(iData),sem(iData)] = weightedStats(currentData(:,1),currentData(:,2),'w');
sd(iData) = sem(iData) * sqrt(sum(currentData(:,2)));
% weighted median: where we're at middle weight
% may need some tweaking
goodCurrentData = sortrows(currentData(all(isfinite(currentData),2),:),1);
weightList = cumsum(goodCurrentData(:,2));
weightList = weightList / weightList(end);
md(iData) = goodCurrentData(find(weightList>0.5,1,'first'),1);
else
m(iData) = nanmean(currentData);
md(iData) = nanmedian(currentData);
sd(iData) = nanstd(currentData);
sem(iData) = sd(iData)/sqrt(sum(isfinite(currentData)));
end

if opt.showMM == 6
% read quantiles - "y"-value, plus x-start-stop
% re-use md array which allows using a loop below instead of
% lots of copy-paste
% md array is md/q1/q3, with third dimension y/xmin/xmax

md(iData,2,1) = prctile(currentData,25);
md(iData,3,1) = prctile(currentData,75);

for qq = 1:3
% find corresponding y-bin
yLoc =  repmat(...
any(yArray>md(iData,qq,1),1) & any(yArray<=md(iData,qq,1),1),...
[4 1]);
% look up corresponding x-values. Note that there is a bit
% of a risk that the line will be exactly between two very
% different bins - but if we make the line longer, it will
% be ugly almost all the time
md(iData,qq,2) = min( xArray( yLoc ) );
md(iData,qq,3) = max( xArray( yLoc ) );
end

end
end
end % loop

sh = [];
if isHistogram
disp('Option addSpread is unavailable if data is supplied as histograms. Call plotSpread separately')
else
try
set(sh{1},'color',[0,128,255]/255);
catch me
if strcmp(me.identifier,'MATLAB:UndefinedFunction')
else
rethrow(me)
end
end
end
end

mh = [];mdh=[];
if opt.showMM
% plot mean, median. Mean is filled red circle, median is green square
% I don't know of a very clever way to flip xy and keep everything
% readable, thus it'll be copy-paste
switch opt.xyOri
case 'normal'
if any(opt.showMM==[1,2])
mh = plot(ah,opt.xValues+xOffset,m,'+r','Color','r','MarkerSize',12);
end
if any(opt.showMM==[1,3])
mdh = plot(ah,opt.xValues+xOffset,md,'sg','MarkerSize',12);
end
if opt.showMM == 4
mh = plot(ah,opt.xValues+xOffset,m,'+r','Color','r','MarkerSize',12);
mdh = myErrorbar(ah,opt.xValues+xOffset,m,sem);
end
if opt.showMM == 5
mh = plot(ah,opt.xValues+xOffset,m,'+r','Color','r','MarkerSize',12);
mdh = myErrorbar(ah,opt.xValues+xOffset,m,sd);
end
if opt.showMM == 6
mdh(1,:) = plot(ah,squeeze(md(:,1,2:3))',repmat(md(:,1,1)',2,1),'color','r','lineWidth',2);%,'lineStyle','--');
mdh(2,:) = plot(ah,squeeze(md(:,2,2:3))',repmat(md(:,2,1)',2,1),'color','r','lineWidth',1);%,'lineStyle','--');
mdh(3,:) = plot(ah,squeeze(md(:,3,2:3))',repmat(md(:,3,1)',2,1),'color','r','lineWidth',1);%,'lineStyle','--');
end
case 'flipped'
if any(opt.showMM==[1,2])
mh = plot(ah,m,opt.xValues+xOffset,'+r','Color','r','MarkerSize',12);
end
if any(opt.showMM==[1,3])
mdh = plot(ah,md,opt.xValues+xOffset,'sg','MarkerSize',12);
end
if opt.showMM == 4
mh = plot(ah,m,opt.xValues+xOffset,'+r','Color','r','MarkerSize',12);
mdh = myErrorbar(ah,m,opt.xValues+xOffset,[sem,NaN(size(sem))]);
end
if opt.showMM == 5
mh = plot(ah,m,opt.xValues+xOffset,'+r','Color','r','MarkerSize',12);
mdh = myErrorbar(ah,m,opt.xValues+xOffset,[sd,NaN(size(sd))]);
end
if opt.showMM == 6
mdh(1,:) = plot(ah,repmat(md(:,1,1)',2,1),squeeze(md(:,1,2:3))','color','r','lineWidth',2);%,'lineStyle','--');
mdh(2,:) = plot(ah,repmat(md(:,2,1)',2,1),squeeze(md(:,2,2:3))','color','r','lineWidth',1);%,'lineStyle','--');
mdh(3,:) = plot(ah,repmat(md(:,3,1)',2,1),squeeze(md(:,3,2:3))','color','r','lineWidth',1);%,'lineStyle','--');
end
end
end

% find extents of x-axis (or y-axis, if flipped)
minX = min(opt.xValues)-stdWidth;
maxX = max(opt.xValues)+stdWidth;

if ~isnan(xAxLim(1))
% we have previous limits
switch opt.xyOri
case 'normal'
minX = min(minX,xAxLim(1));
maxX = max(maxX,xAxLim(2));
case 'flipped'
minX = min(minX,yAxLim(1));
maxX = max(maxX,yAxLim(2));
end
end

% if ~empty, use xNames
switch opt.xyOri
case 'normal'
switch opt.xMode
case 'manual'
if newAx == false
warning('DISTRIBUTIONPLOT:ERASINGLABELS','Plotting into an existing axes and specifying labels will erase previous labels')
end
set(ah,'XTick',opt.xValues);
if ~isempty(opt.xNames)
set(ah,'XTickLabel',opt.xNames)
end
case 'auto'
% no need to do anything
end
if ~isempty(opt.yLabel)
ylabel(ah,opt.yLabel);
end
% have plot start/end properly
xlim([minX,maxX])
case 'flipped'
switch opt.xMode
case 'manual'
if newAx == false
warning('DISTRIBUTIONPLOT:ERASINGLABELS','Plotting into an existing axes and specifying labels will erase previous labels')
end
set(ah,'YTick',opt.xValues);
if ~isempty(opt.xNames)
set(ah,'YTickLabel',opt.xNames)
end
case 'auto'
% no need to do anything
end
if ~isempty(opt.yLabel)
xlabel(ah,opt.yLabel);
end
% have plot start/end properly
ylim([minX,maxX])
end

%==========================

%==========================
%% CLEANUP & ASSIGN OUTPUT
%==========================

if nargout > 0
handles{1} = hh;
handles{2} = [mh;mdh];
handles{3} = ah;
handles{4} = sh;
end

set(ah,'NextPlot',holdState);```