Code covered by the BSD License  

Highlights from
Violin Plots for plotting multiple distributions (distributionPlot.m)

image thumbnail

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.
%                   Please download plotSpread.m separately from the File
%                   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)
%         distributionPlot({data(1:5:end),repmat(data,2,1)},'addSpread',true,'showMM',false,'histOpt',2) 
%
%         %--- 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
%           hh=distributionPlot({xx,xx,xx},'color','g','addSpread',1,'histOpt',2,'showMM',0);
%           set(hh{4}{1},'color','r','marker','o')
% END
%
% REMARKS To show distributions as clouds of points (~beeswarm plot),
%         and/or to use the option "addSpread", please download the
%         additional function plotSpread.m from the File Exchange
%         http://www.mathworks.com/matlabcentral/fileexchange/37105-plot-spread-points-beeswarm-plot
%         
%         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.
%           handles = distributionPlot(data,distWidth,showMM,xNames,histOpt,divFactor,invert,addSpread,globalNorm)
%           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
%
%
% SEE ALSO histogram, ksdensity, plotSpread, boxplot, grp2idx
%

% 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.addSpread = false;
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 any(badCol)
        nCols = cellfun(@(x)(size(x,2)),data(badCol));
        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',...
                num2str(find(badCol)'));
            data(badCol) = cellfun(@(x)(x(:)),data(badCol),'UniformOutput',false);
        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
    % distWidth,showMM,xNames,histOpt,divFactor,invert,addSpread,globalNorm
    def.distWidth = 1.9;
    parserObj.addOptional('distWidth',def.distWidth);
    parserObj.addOptional('showMM',def.showMM);
    parserObj.addOptional('xNames',def.xNames);
    parserObj.addOptional('histOpt',def.histOpt);
    parserObj.addOptional('divFactor',def.divFactor);
    parserObj.addOptional('invert',def.invert);
    parserObj.addOptional('addSpread',def.addSpread);
    parserObj.addOptional('globalNorm',def.globalNorm);
    parserObj.addOptional('groups',def.groups);
    parserObj.addOptional('yLabel',def.yLabel);
    parserObj.addOptional('color',def.color);
    
    
    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(:)'
        parserObj.addParamValue(dn{1},def.(dn{1}));
    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');
set(ah,'NextPlot','add');

% 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;
    
    % add patch
    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 opt.addSpread
    if isHistogram
        disp('Option addSpread is unavailable if data is supplied as histograms. Call plotSpread separately')
    else
        % add spread
        try
            sh = plotSpread(ah,data,'xValues',opt.xValues,'xyOri',opt.xyOri);
            set(sh{1},'color',[0,128,255]/255);
        catch me
            if strcmp(me.identifier,'MATLAB:UndefinedFunction')
                error('plotSpread not found. Please download it from the Matlab File Exchange')
            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);

Contact us