--- boxplot_v2.m	2020-06-06 00:27:32.205253700 +0200
+++ boxplot_v3.m	2020-06-06 00:24:56.604803600 +0200
@@ -25,100 +25,117 @@
 ## @deftypefnx {Function File} {@var{s} =} boxplot (@var{data}, @var{group}, @var{notched}, @var{symbol}, @var{orientation}, @var{whisker}, @dots{})
 ## @deftypefnx {Function File} {@var{s} =} boxplot (@var{data}, @var{options})
 ## @deftypefnx {Function File} {@var{s} =} boxplot (@var{data}, @var{group}, @var{options}, @dots{})
-##
-## @deftypefnx {Function File} {[@dots{} @var{h}]=} boxplot (@dots{})
-##
+## @deftypefnx {Function File} {[@dots{} @var{h}]=} boxplot (@var{data}, @dots{})
 ## Produce a box plot.
 ##
-## The box plot is a graphical display that simultaneously describes several
+## A box plot is a graphical display that simultaneously describes several
 ## important features of a data set, such as center, spread, departure from
 ## symmetry, and identification of observations that lie unusually far from
 ## the bulk of the data.
 ##
+## Input arguments recognized by boxplot are:
+##
+## @itemize
+## @item
 ## @var{data} is a matrix with one column for each data set, or a cell vector
 ## with one cell for each data set.  Each cell must contain a numerical row or
 ## column vector (NaN and NA are ignored) and not a nested vector of cells.
 ##
+## @item
 ## @var{notched} = 1 produces a notched-box plot.  Notches represent a robust
 ## estimate of the uncertainty about the median.
 ##
 ## @var{notched} = 0 (default) produces a rectangular box plot.
 ##
-## @var{notched} in (0,1) produces a notch of the specified depth.
-## notched values outside (0,1) are amusing if not exactly practical.
+## @var{notched} within the interval (0,1) produces a notch of the specified
+## depth.  Notched values outside (0,1) are amusing if not exactly impractical.
 ##
-## @var{symbol} sets the symbol for the outlier values, default symbol for
-## points that lie outside 3 times the interquartile range is 'o',
-## default symbol for points between 1.5 and 3 times the interquartile
-## range is '+'.
+## @item
+## @var{symbol} sets the symbol for the outlier values.  The default symbol
+## for points that lie outside 3 times the interquartile range is 'o';
+## the default symbol for points between 1.5 and 3 times the interquartile
+## range is '+'. @*
+## Alternative settings:
 ##
-## @var{symbol} = '.' points between 1.5 and 3 times the IQR is marked with
+## @var{symbol} = '.': points between 1.5 and 3 times the IQR are marked with
 ## '.' and points outside 3 times IQR with 'o'.
 ##
-## @var{symbol} = ['x','*'] points between 1.5 and 3 times the IQR is marked with
-## 'x' and points outside 3 times IQR with '*'.
-##
-## @var{orientation} = 0 makes the boxes horizontal, by default @var{orientation} = 1,
-## which plots the boxes vertically.  Alternatively, options can be  passed as a string.
+## @var{symbol} = ['x','*']: points between 1.5 and 3 times the IQR are marked
+## with 'x' and points outside 3 times IQR with '*'.
 ##
-## @var{orientation} = 'vertical'      (default value)
-##
-## @var{orientation} = 'horizontal'
+## @item
+## @var{orientation} = 0 makes the boxes horizontal (default).
+## @var{orientation} = 1 plots the boxes vertically.  Alternatively,
+## orientation options can be passed as a string, e.g., 'vertical' or
+## 'horizontal'.
 ##
+## @item
 ## @var{whisker} defines the length of the whiskers as a function of the IQR
-## (default = 1.5). If @var{whisker} = 0 then @code{boxplot} displays all data
+## (default = 1.5).  If @var{whisker} = 0 then @code{boxplot} displays all data
 ## values outside the box using the plotting symbol for points that lie
 ## outside 3 times the IQR.
 ##
-## @var{group} may be passed as an optional argument only in the second position after
-## @var{data} containing a numerical vector, which defines separate categories, each
-## plotted in a different box, for each set of @var{DATA} values that share the same
-## @var{group} value or values.  With the formalism (@var{data}, @var{group}), both must
-## be vectors of the same length.
-##
-## @var{options} are additional paired arguments passed with the formalism (Name, Value)
-## that provide extra functionality as listed below.  @var{options} can be passed at any
-## order after the initial arguments.
+## @item
+## @var{group} may be passed as an optional argument only in the second
+## position after @var{data}.  @var{group} contains a numerical vector defining
+## separate categories, each plotted in a different box, for each set of
+## @var{DATA} values that share the same @var{group} value or values.  With
+## the formalism (@var{data}, @var{group}), both must be vectors of the same
+## length.
+##
+## @item
+## @var{options} are additional paired arguments passed with the formalism
+## (Name, Value) that provide extra functionality as listed below.
+## @var{options} can be passed at any order after the initial arguments.
 ##
 ## @multitable {Name} {Value} {description} @columnfractions .2 .2 .6
-## @item 'Notch' @tab  'on' @tab notched by 0.25 of the boxes width
-## @item @tab 'off' @tab produces a straight box
-## @item @tab scalar @tab proportional width of the notch
+## @item 'Notch' @tab  'on' @tab Notched by 0.25 of the boxes width.
+## @item @tab 'off' @tab Produces a straight box.
+## @item @tab scalar @tab Proportional width of the notch.
 ##
-## @item 'Symbol' @tab '.' @tab defines only outliers between 1.5 and 3 IQR
+## @item 'Symbol' @tab '.' @tab Defines only outliers between 1.5 and 3 IQR.
 ## @item @tab ['x','*'] @tab 2nd character defines outliers > 3 IQR
 ##
-## @item 'Orientation' @tab 'vertical' @tab default value, can also be defined with 1
-## @item @tab 'horizontal' @tab can also be defined with 0
-##
-## @item 'Whisker' @tab scalar @tab multiplier of IQR (defualt is 1.5)
-##
-## @item 'OutlierTags' @tab 'on' @tab plot the vector index of the outlier value next to its point
-## @item @tab 'off' @tab no tags are plotted (default value)
-##
-## @item 'Sample_IDs' @tab 'cell' @tab a cell vector with one cell for each data set containing
-## a nested cell vector with each sample's ID (should be a string).  If this option is passed, then
-## all outliers are tagged with their respective sample's id string instead of their vector's index
-##
-## @item 'BoxWidth' @tab 'proportional' @tab create boxes with their width proportional to the number
-## of samples in their respective dataset (default value)
-## @item @tab 'fixed' @tab make all boxes with equal width
-##
-## @item 'Widths' @tab scalar @tab scaling factor for box widths (default value is 0.4)
-##
-## @item 'BoxStyle' @tab 'outline' @tab draw boxes as outlines (default value)
-## @item @tab 'filled' @tab fill boxes with a color (outlines are still plotted)
-##
-## @item 'Positions' @tab vector @tab numerical vector that defines the position of each data set.
-## It must have the same length as the number of groups and can be used to reorder or group them
-## in a desired manner.  This vector merely defines the points along the group axis, which by default
-## is [1:number of groups].
-##
-## @item 'Labels' @tab cell @tab a cell vector of strings containing the names of each group.  By
-## default each group is labeled numerically according to its order in the data set
+## @item 'Orientation' @tab 'vertical' @tab Default value, can also be defined
+## with numerical 1.
+## @item @tab 'horizontal' @tab Can also be defined with numerical 0.
+##
+## @item 'Whisker' @tab scalar @tab Multiplier of IQR (defualt is 1.5).
+##
+## @item 'OutlierTags' @tab 'on' @tab Plot the vector index of the outlier
+## value next to its point.
+## @item @tab 'off' @tab No tags are plotted (default value).
+##
+## @item 'Sample_IDs' @tab 'cell' @tab A cell vector with one cell for each
+## data set containing a nested cell vector with each sample's ID (should be
+## a string).  If this option is passed, then all outliers are tagged with
+## their respective sample's ID string instead of their vector's index.
+##
+## @item 'BoxWidth' @tab 'proportional' @tab Create boxes with their width
+## proportional to the number of samples in their respective dataset (default
+## value).
+## @item @tab 'fixed' @tab Make all boxes with equal width.
+##
+## @item 'Widths' @tab scalar @tab Scaling factor for box widths (default
+## value is 0.4).
+##
+## @item 'BoxStyle' @tab 'outline' @tab Draw boxes as outlines (default value).
+## @item @tab 'filled' @tab Fill boxes with a color (outlines are still
+## plotted).
+##
+## @item 'Positions' @tab vector @tab Numerical vector that defines the
+## position of each data set.  It must have the same length as the number of
+## in a desired manner.  This vector merely defines the points along the group
+## axis, which by default is [1:number of groups].
+##
+## @item 'Labels' @tab cell @tab A cell vector of strings containing the names
+## of each group.  By default each group is labeled numerically according to
+## its order in the data set
 ## @end multitable
+## @end itemize
 ##
-## Supplemental arguments (@dots{}) are concatenated and passed to plot.
+## Supplemental arguments not described above (@dots{}) are concatenated and
+## passed to the plot() function.
 ##
 ## The returned matrix @var{s} has one column for each data set as follows:
 ##
@@ -132,8 +149,8 @@
 ## @item 7 @tab Upper confidence limit for median
 ## @end multitable
 ##
-## The returned structure @var{h} has handles to the plot elements, allowing
-## customization of the visualization using set/get functions.
+## The returned structure @var{h} contains handles to the plot elements,
+## allowing customization of the visualization using set/get functions.
 ##
 ## Example
 ##
@@ -146,46 +163,47 @@
 ##
 ## @end deftypefn
 
-function [s hs] = boxplot (data, varargin)
+function [s_o, hs_o] = boxplot (data, varargin)
 
-  ## assign parameter defaults
+  ## Assign parameter defaults
   if (nargin < 1)
     print_usage;
   endif
 
-  ## default values
+  ## Default values
   maxwhisker = 1.5;
   orientation = 1;
-  symbol = ['+', 'o'];
+  symbol = ["+", "o"];
   notched = 0;
   plot_opts = {};
   groups = [];
   sample_IDs = {};
   outlier_tags = 0;
-  box_width = 'proportional';
+  box_width = "proportional";
   widths = 0.4;
-  box_style = 0;	
+  box_style = 0;
   positions = [];
   labels = {};
 
   ## Optional arguments analysis
   numarg = nargin - 1;
-  option_args = {'Notch'; 'Symbol'; 'Orientation'; 'Whisker'; 'OutlierTags';
-                'Sample_IDs'; 'BoxWidth'; 'Widths'; 'BoxStyle'; 'Positions'; 'Labels'};
+  option_args = {"Notch"; "Symbol"; "Orientation"; "Whisker"; ...
+                 "OutlierTags"; "Sample_IDs"; "BoxWidth"; "Widths"; ...
+                 "BoxStyle"; "Positions"; "Labels"};
   indopt = 1;
   group_exists = 0;
   while (numarg)
     dummy = varargin{indopt++};
     if (! ischar (dummy) && indopt < 6)
-      ## MatLAB allows passing the second argument as a grouping vector
+      ## MATLAB allows passing the second argument as a grouping vector
       if (length (dummy) > 1)
         if (2 != indopt)
-          error ('Boxplot.m: grouping vector may only be passed as second arg');
+          error ("Boxplot: grouping vector may only be passed as second arg");
         endif
         groups = dummy;
         group_exists = 1;
       elseif (length (dummy) == 1)
-        ## old way: positional argument
+        ## Old way: positional argument
         switch indopt - group_exists
           case 2
             notched = dummy;
@@ -200,43 +218,43 @@
       numarg--;
       continue;
     else
-      if (3 == indopt - group_exists && length (dummy) <= 2)
+      if (3 == (indopt - group_exists) && length (dummy) <= 2)
         symbol = dummy;
         numarg--;
         continue;
       else
-        ## checking for additional paired arguments
-        tt = find(strcmp(dummy, option_args));
+        ## Check for additional paired arguments
+        tt = find (strcmpi (dummy, option_args));
         switch (tt)
           case 1
             notched = varargin{indopt};
-            ## chech for string input: 'on' or 'off'
-            if strcmp (notched, "on")
+            ## Check for string input: "on" or "off"
+            if strcmpi (notched, "on")
               notched = 1;
             endif
-            if strcmp (notched, "off")
+            if strcmpi (notched, "off")
               notched = 0;
             endif
           case 2
             symbol = varargin{indopt};
           case 3
             orientation = varargin{indopt};
-            ## chech for string input: 'vertical' or 'horizontal'
-            if strcmp (notched, "vertical")
+            ## Check for string input: "vertical" or "horizontal"
+            if strcmpi (notched, "vertical")
               notched = 1;
             endif
-            if strcmp (notched, "horizontal")
+            if strcmpi (notched, "horizontal")
               notched = 0;
             endif
           case 4
             maxwhisker = varargin{indopt};
           case 5
             outlier_tags = varargin{indopt};
-            ## chech for string input: 'on' or 'off'
-            if strcmp (outlier_tags, "on")
+            ## Check for string input: "on" or "off"
+            if strcmpi (outlier_tags, "on")
               outlier_tags = 1;
             endif
-            if strcmp (outlier_tags, "off")
+            if strcmpi (outlier_tags, "off")
               outlier_tags = 0;
             endif
           case 6
@@ -244,21 +262,22 @@
             outlier_tags = 1;
           case 7
             box_width = varargin{indopt};
-            ## check for string input: 'fixed' or 'proportional' (default if mis-spelled)
-            if strcmp (box_width, "fixed")
-              box_width = 'fixed';
+            ## Check for string input: "fixed" or "proportional" (default
+            ## if mis-spelled)
+            if strcmpi (box_width, "fixed")
+              box_width = "fixed";
             else
-              box_width = 'proportional';
+              box_width = "proportional";
             endif
           case 8
             widths = varargin{indopt};
           case 9
             box_style = varargin{indopt};
-            ## chech for string input: 'outline' or 'filled'
-            if strcmp (box_style, "outline")
+            ## Check for string input: "outline" or "filled"
+            if strcmpi (box_style, "outline")
               box_style = 0;
             endif
-            if strcmp (box_style, "filled")
+            if strcmpi (box_style, "filled")
               box_style = 1;
             endif
           case 10
@@ -266,12 +285,12 @@
           case 11
             labels = varargin{indopt};
           otherwise
-            ## take two args and append them to plot_opts
-            plot_opts(1, end+1:end+2) = {dummy,  varargin{indopt}};
+            ## Take two args and append them to plot_opts
+            plot_opts(1, end+1:end+2) = {dummy, varargin{indopt}};
         endswitch
       endif
       numarg -= 2;
-      indopt++; 
+      indopt++;
     endif
   endwhile
 
@@ -284,14 +303,14 @@
   endif
   a = 1-notched;
 
-  ## figure out how many data sets we have
+  ## Figure out how many data sets we have
   if (isempty (groups))
     if (iscell (data))
       nc = length (data);
       for ind_c = (1:nc)
         lc(ind_c) = length (data{ind_c});
       endfor
-	  else
+    else
       if (isvector (data))
         data = data(:);
       endif
@@ -299,30 +318,30 @@
       lc = ones (1,nc) * rows (data);
     endif
     groups = (1:nc);
-    ## check if sample_IDs exists that it has same size with data
+    ## Check if sample_IDs exists that it has same size with data
     if (! isempty (sample_IDs) && length (sample_IDs) == 1)
       for ind_c = (1:nc)
         if (lc(ind_c) != length (sample_IDs))
-          error ('Boxplot.m: Sample_IDs must match the data');
+          error ("Boxplot: Sample_IDs must match the data");
         endif
       endfor
     elseif (! isempty (sample_IDs) && length (sample_IDs) == nc)
       for ind_c = (1:nc)
         if (lc(ind_c) != length (sample_IDs{ind_c}))
-          error ('Boxplot.m: Sample_IDs must match the data');
+          error ("Boxplot: Sample_IDs must match the data");
         endif
       endfor
     elseif (! isempty (sample_IDs) && length (sample_IDs) != nc)
-      error ('Boxplot.m: Sample_IDs must match the data');
+      error ("Boxplot: Sample_IDs must match the data");
     endif
   else
     if (! isvector (data))
-      error ('Boxplot.m: with the formalism (data, group), both must be vectors');
+      error ("Boxplot: with the formalism (data, group), both must be vectors");
     end
-    ## check if sample IDs exist that they have same size with data
+    ## Check if sample IDs exist that they have same size with data
     if (! isempty (sample_IDs))
       if (length (sample_IDs) != 1 || length (sample_IDs{1}) != length (data))
-        error ('Boxplot.m: Sample_IDs must match the data');
+        error ("Boxplot: Sample_IDs must match the data");
       endif
       nc = unique (groups);
       dummy_data = cell (1, length (nc));
@@ -346,14 +365,14 @@
     endif
   endif
 
-  ## compute statistics
+  ## Compute statistics.
   ## s will contain
   ##    1,5    min and max
   ##    2,3,4  1st, 2nd and 3rd quartile
   ##    6,7    lower and upper confidence intervals for median
   s = zeros (7, nc);
   box = zeros (1, nc);
-  ## create labels according to number of datasets as ordered in parsed data
+  ## Create labels according to number of datasets as ordered in parsed data
   ## in case they are not provided by the user as optional argument
   if (isempty (labels))
     for i=1:nc
@@ -361,13 +380,15 @@
       labels(i) = {column_label};
     endfor
   endif
-  ## arrange the boxes into desired positions (if requested, otherwise leave default 1:nc)
+  ## Arrange the boxes into desired positions (if requested, otherwise leave
+  ## default 1:nc)
   if (! isempty (positions))
     groups = positions;
-  endif	
-  ## initialize whisker matrices to correct size and all necessary outlier variables
-  whisker_x = ones (2,1)*[groups, groups];
-  whisker_y = zeros (2, 2*nc);
+  endif
+  ## Initialize whisker matrices to correct size and all necessary outlier
+  ## variables
+  whisker_x = ones (2, 1) * [groups, groups];
+  whisker_y = zeros (2, 2 * nc);
   outliers_x = [];
   outliers_y = [];
   outliers_idx = [];
@@ -390,8 +411,8 @@
       col = data(:, indi);
       sIDs = num2cell([1:length(col)]);
     endif
-    ## Skip missing data (NaN, NA) and remove respective sample IDs
-    ## do this only on nonempty data
+    ## Skip missing data (NaN, NA) and remove respective sample IDs.
+    ## Do this only on nonempty data
     if length(col) > 0
       remove_samples = find (col(isnan (col) | isna (col)));
       if length (remove_samples) > 0
@@ -403,88 +424,90 @@
     nd = length (col);
     box(indi) = nd;
     if (nd > 1)
-      ## min,max and quartiles
+      ## Min, max and quartiles
       s(1:5, indi) = statistics (col)(1:5);
-      ## confidence interval for the median
-      est = 1.57*(s(4, indi)-s(2, indi))/sqrt (nd);
-      s(6, indi) = max ([s(3, indi)-est, s(2, indi)]);
-      s(7, indi) = min ([s(3, indi)+est, s(4, indi)]);
-      ## whiskers out to the last point within the desired inter-quartile range
-      IQR = maxwhisker*(s(4, indi)-s(2, indi));
-      whisker_y(:, indi) = [min(col(col >= s(2, indi)-IQR)); s(2, indi)];
-      whisker_y(:,nc+indi) = [max(col(col <= s(4, indi)+IQR)); s(4, indi)];
-      ## outliers beyond 1 and 2 inter-quartile ranges
-      outliers = col((col < s(2, indi)-IQR & col >= s(2, indi)-2*IQR) | (col > s(4, indi)+IQR & col <= s(4, indi)+2*IQR));
-      outliers2 = col(col < s(2, indi)-2*IQR | col > s(4, indi)+2*IQR);
-      ## get outliers indices from this dataset
-      if length (outliers) > 0
+      ## Confidence interval for the median
+      est = 1.57 * (s(4, indi) - s(2, indi)) / sqrt (nd);
+      s(6, indi) = max ([s(3, indi) - est, s(2, indi)]);
+      s(7, indi) = min ([s(3, indi) + est, s(4, indi)]);
+      ## Whiskers out to the last point within the desired inter-quartile range
+      IQR = maxwhisker * (s(4, indi) - s(2, indi));
+      whisker_y(:, indi) = [min(col(col >= s(2, indi) - IQR)); s(2, indi)];
+      whisker_y(:, nc+indi) = [max(col(col <= s(4, indi) + IQR)); s(4, indi)];
+      ## Outliers beyond 1 and 2 inter-quartile ranges
+      outliers = col((col < s(2, indi) - IQR & col >= s(2, indi) - 2 * IQR) | ...
+                     (col > s(4, indi) + IQR & col <= s(4, indi) + 2 * IQR));
+      outliers2 = col(col < s(2, indi) - 2 * IQR | col > s(4, indi) + 2 * IQR);
+      ## Get outliers indices from this dataset
+      if (length (outliers) > 0)
         for out_i = 1:length (outliers)
-          outliers_idx = [outliers_idx; find(col == outliers(out_i))];
-          outliers_IDs = {outliers_IDs{}, sIDs{find(col == outliers(out_i))}};
+          outliers_idx = [outliers_idx; (find (col == outliers(out_i)))];
+          outliers_IDs = {outliers_IDs{:}, sIDs{(find (col == outliers(out_i)))}};
         endfor
       endif
-      if length (outliers2) > 0
+      if (length (outliers2) > 0)
         for out_i = 1:length (outliers2)
           outliers2_idx = [outliers2_idx; find(col == outliers2(out_i))];
-          outliers2_IDs = {outliers2_IDs{}, sIDs{find(col == outliers2(out_i))}};
+          outliers2_IDs = {outliers2_IDs{:}, sIDs{find(col == outliers2(out_i))}};
         endfor
       endif
-      outliers_x = [outliers_x; groups(indi)*ones(size(outliers))];
+      outliers_x = [outliers_x; (groups(indi) * ones (size (outliers)))];
       outliers_y = [outliers_y; outliers];
-      outliers2_x = [outliers2_x; groups(indi)*ones(size(outliers2))];
+      outliers2_x = [outliers2_x; (groups(indi) * ones (size (outliers2)))];
       outliers2_y = [outliers2_y; outliers2];
     elseif (1 == nd)
-      ## all statistics collapse to the value of the point
+      ## All statistics collapse to the value of the point
       s(:, indi) = col;
-      ## single point data sets are plotted as outliers.
+      ## Single point data sets are plotted as outliers.
       outliers_x = [outliers_x; groups(indi)];
       outliers_y = [outliers_y; col];
-      ## append the single point's index to keep the outliers' vector aligned
+      ## Append the single point's index to keep the outliers' vector aligned
       outliers_idx = [outliers_idx; 1];
-      outliers_IDs = {outliers_IDs{}, sIDs{}};
+      outliers_IDs = {outliers_IDs{:}, sIDs{:}};
     else
-      ## no statistics if no points
+      ## No statistics if no points
       s(:, indi) = NaN;
-    end
-  end
+    endif
+  endfor
 
   ## Note which boxes don't have enough stats
   chop = find (box <= 1);
 
-  ## replicate widths (if scalar or shorter vector) to match the number of boxes
-  widths = widths(repmat(1:length(widths),1,nc));
-  ## truncate just in case :)
+  ## Replicate widths (if scalar or shorter vector) to match the number of boxes
+  widths = widths(repmat (1:length (widths), 1, nc));
+  ## Truncate just in case :)
   widths([nc+1:end]) = [];
   ## Draw a box around the quartiles, with box width being fixed or proportional
   ## to the number of items in the box.
-  if (strcmp (box_width, 'proportional'))
+  if (strcmpi (box_width, "proportional"))
     box = box .* (widths ./ max(box));
   else
     box = box .* (widths ./ box);
   endif
   ## Draw notches if desired.
-  quartile_x = ones (11,1)*groups + [-a;-1;-1;1;1;a;1;1;-1;-1;-a]*box;
-  quartile_y = s([3,7,4,4,7,3,6,2,2,6,3],:);
+  quartile_x = ones (11, 1) * groups + ...
+               [-a; -1; -1; 1 ; 1; a; 1; 1; -1; -1; -a] * box;
+  quartile_y = s([3, 7, 4, 4, 7, 3, 6, 2, 2, 6, 3], :);
 
   ## Draw a line through the median
-  median_x = ones (2,1)*groups + [-a;+a]*box;
-  median_y = s([3,3],:);
+  median_x = ones (2, 1) * groups + [-a; +a] * box;
+  median_y = s([3, 3], :);
 
   ## Chop all boxes which don't have enough stats
   quartile_x(:, chop) = [];
   quartile_y(:, chop) = [];
-  whisker_x(:,[chop, chop+nc]) = [];
-  whisker_y(:,[chop, chop+nc]) = [];
+  whisker_x(:, [chop, chop + nc]) = [];
+  whisker_y(:, [chop, chop + nc]) = [];
   median_x(:, chop) = [];
   median_y(:, chop) = [];
 
   ## Add caps to the remaining whiskers
   cap_x = whisker_x;
-  cap_x(1, :) -= 0.05;
-  cap_x(2, :) += 0.05;
+  cap_x(1, :) -= repmat ((widths / 1.5), 1, 2);
+  cap_x(2, :) += repmat ((widths / 1.5), 1, 2);
   cap_y = whisker_y([1, 1], :);
 
-  ## calculate coordinates for outlier tags
+  ## Calculate coordinates for outlier tags
   outliers_tags_x = outliers_x + 0.08;
   outliers_tags_y = outliers_y;
   outliers2_tags_x = outliers2_x + 0.08;
@@ -492,7 +515,7 @@
 
   ## Do the plot
   if (orientation)
-    ## define outlier_tags' vertical alignment
+    ## Define outlier_tags' vertical alignment
     outlier_tags_alignment = {"horizontalalignment", "left"};
     if (isempty (plot_opts))
       if (box_style)
@@ -500,19 +523,19 @@
         hold on;
       endif
       h = plot (quartile_x, quartile_y, "b;;",
-          whisker_x, whisker_y, "b;;",
-          cap_x, cap_y, "b;;",
-          median_x, median_y, "r;;",
-          outliers_x, outliers_y, [symbol(1), "r;;"],
-          outliers2_x, outliers2_y, [symbol(2), "r;;"]);
-      ## print outlier tags
+                whisker_x, whisker_y, "b;;",
+                cap_x, cap_y, "b;;",
+                median_x, median_y, "r;;",
+                outliers_x, outliers_y, [symbol(1), "r;;"],
+                outliers2_x, outliers2_y, [symbol(2), "r;;"]);
+      ## Print outlier tags
       if (outlier_tags == 1 && outliers_x > 0)
         t1 = plot_tags (outliers_tags_x, outliers_tags_y, outliers_idx,
-             outliers_IDs, sample_IDs, outlier_tags_alignment);
+                        outliers_IDs, sample_IDs, outlier_tags_alignment);
       endif
       if (outlier_tags == 1 && outliers2_x > 0)
         t2 = plot_tags (outliers2_tags_x, outliers2_tags_y, outliers2_idx,
-             outliers2_IDs, sample_IDs, outlier_tags_alignment);
+                        outliers2_IDs, sample_IDs, outlier_tags_alignment);
       endif
     else
       if (box_style)
@@ -520,88 +543,96 @@
         hold on;
       endif
       h = plot (quartile_x, quartile_y, "b;;",
-          whisker_x, whisker_y, "b;;",
-          cap_x, cap_y, "b;;",
-          median_x, median_y, "r;;",
-          outliers_x, outliers_y, [symbol(1), "r;;"],
-          outliers2_x, outliers2_y, [symbol(2), "r;;"], plot_opts{:});
-      ## print outlier tags
+                whisker_x, whisker_y, "b;;",
+                cap_x, cap_y, "b;;",
+                median_x, median_y, "r;;",
+                outliers_x, outliers_y, [symbol(1), "r;;"],
+                outliers2_x, outliers2_y, [symbol(2), "r;;"], plot_opts{:});
+      ## Print outlier tags
       if (outlier_tags == 1 && outliers_x > 0)
         t1 = plot_tags (outliers_tags_x, outliers_tags_y, outliers_idx,
-             outliers_IDs, sample_IDs, outlier_tags_alignment);
+                        outliers_IDs, sample_IDs, outlier_tags_alignment);
       endif
       if (outlier_tags == 1 && outliers2_x > 0)
         t2 = plot_tags (outliers2_tags_x, outliers2_tags_y, outliers2_idx,
-             outliers2_IDs, sample_IDs, outlier_tags_alignment);
+                        outliers2_IDs, sample_IDs, outlier_tags_alignment);
       endif
     endif
   else
-    ## define outlier_tags' horizontal alignment
+   ## Define outlier_tags' horizontal alignment
     outlier_tags_alignment = {"horizontalalignment", "left", "rotation", 90};
     if (isempty (plot_opts))
+      h = plot (quartile_y, quartile_x, "b;;",
+               whisker_y, whisker_x, "b;;",
+               cap_y, cap_x, "b;;",
+               median_y, median_x, "r;;",
+               outliers_y, outliers_x, [symbol(1), "r;;"],
+               outliers2_y, outliers2_x, [symbol(2), "r;;"]);
       if (box_style)
-        f = fill(quartile_y, quartile_x,"y");
+        f = fill (quartile_y, quartile_x, "y");
         hold on;
       endif
       h = plot (quartile_y, quartile_x, "b;;",
-          whisker_y, whisker_x, "b;;",
-          cap_y, cap_x, "b;;",
-          median_y, median_x, "r;;",
-          outliers_y, outliers_x, [symbol(1), "r;;"],
-          outliers2_y, outliers2_x, [symbol(2), "r;;"]);
-      ## print outlier tags
+                whisker_y, whisker_x, "b;;",
+                cap_y, cap_x, "b;;",
+                median_y, median_x, "r;;",
+                outliers_y, outliers_x, [symbol(1), "r;;"],
+                outliers2_y, outliers2_x, [symbol(2), "r;;"]);
+      ## Print outlier tags
       if (outlier_tags == 1 && outliers_x > 0)
         t1 = plot_tags (outliers_tags_y, outliers_tags_x, outliers_idx,
-             outliers_IDs, sample_IDs, outlier_tags_alignment);
+                        outliers_IDs, sample_IDs, outlier_tags_alignment);
       endif
       if (outlier_tags == 1 && outliers2_x > 0)
         t2 = plot_tags (outliers2_tags_y, outliers2_tags_x, outliers2_idx,
-             outliers2_IDs, sample_IDs, outlier_tags_alignment);
+                        outliers2_IDs, sample_IDs, outlier_tags_alignment);
       endif
     else
       if (box_style)
-        f = fill(quartile_y, quartile_x,"y");
+        f = fill (quartile_y, quartile_x,"y");
         hold on;
       endif
       h = plot (quartile_y, quartile_x, "b;;",
-          whisker_y, whisker_x, "b;;",
-          cap_y, cap_x, "b;;",
-          median_y, median_x, "r;;",
-          outliers_y, outliers_x, [symbol(1), "r;;"],
-          outliers2_y, outliers2_x, [symbol(2), "r;;"], plot_opts{:});
-      ## print outlier tags
+                whisker_y, whisker_x, "b;;",
+                cap_y, cap_x, "b;;",
+                median_y, median_x, "r;;",
+                outliers_y, outliers_x, [symbol(1), "r;;"],
+                outliers2_y, outliers2_x, [symbol(2), "r;;"], plot_opts{:});
+      ## Print outlier tags
       if (outlier_tags == 1 && outliers_x > 0)
         t1 = plot_tags (outliers_tags_y, outliers_tags_x, outliers_idx,
-             outliers_IDs, sample_IDs, outlier_tags_alignment);
+                        outliers_IDs, sample_IDs, outlier_tags_alignment);
       endif
       if (outlier_tags == 1 && outliers2_x > 0)
         t2 = plot_tags (outliers2_tags_y, outliers2_tags_x, outliers2_idx,
-             outliers2_IDs, sample_IDs, outlier_tags_alignment);
+                        outliers2_IDs, sample_IDs, outlier_tags_alignment);
       endif
     endif
   endif
 
   ## Distribute handles for box outlines and box fill (if any)
-  nq = 1:size(quartile_x,2);
+  nq = 1 : size(quartile_x, 2);
   hs.box = h(nq);
   if (box_style)
-    nf = 1:length(groups);
+    nf = 1 : length(groups);
     hs.box_fill = f(nf);
   else
     hs.box_fill = [];
   endif
+
   ## Distribute handles for whiskers (including caps) and median lines
-  nw = nq(end) + [1:2*size(whisker_x,2)];
+  nw = nq(end) + [1 : 2 * (size (whisker_x, 2))];
   hs.whisker = h(nw);
-  nm = nw(end)+ [1:size(median_x,2)];
+  nm = nw(end)+ [1 : (size (median_x, 2))];
   hs.median = h(nm);
-  ## Distribute handles for outliers (if any) and their respective tags (if applicable)
+  ## Distribute handles for outliers (if any) and their respective tags
+  ## (if applicable)
   no = nm;
   if (! isempty (outliers_y))
-    no = nm(end) + [1:size(outliers_y,2)];
+    no = nm(end) + [1 : size(outliers_y, 2)];
     hs.outliers = h(no);
     if (outlier_tags == 1)
-      nt = 1:length(outliers_tags_y);
+      nt = 1 : length (outliers_tags_y);
       hs.out_tags = t1(nt);
     else
       hs.out_tags = [];
@@ -610,12 +641,13 @@
     hs.outliers = [];
     hs.out_tags = [];
   endif
-  ## Distribute handles for extreme outliers (if any) and their respective tags (if applicable)
+  ## Distribute handles for extreme outliers (if any) and their respective tags
+  ## (if applicable)
   if (! isempty (outliers2_y))
-    no2 = no(end) + [1:size(outliers2_y,2)];
+    no2 = no(end) + [1 : size(outliers2_y, 2)];
     hs.outliers2 = h(no2);
     if (outlier_tags == 1)
-      nt2 = 1:length(outliers2_tags_y);
+      nt2 = 1 : length (outliers2_tags_y);
       hs.out_tags2 = t2(nt2);
     else
       hs.out_tags2 = [];
@@ -625,11 +657,12 @@
     hs.out_tags2 = [];
   end
 
-  ## Redraw the median lines to avoid colour overlapping in case of 'filled' BoxStyle
+  ## Redraw the median lines to avoid colour overlapping in case of 'filled'
+  ## BoxStyle
   if (box_style)
-    set (hs.median,"color","r");
+    set (hs.median, "color", "r");
   endif
-  
+
   ## Print labels according to orientation and return handle
   if (orientation)
     set (gca(), "xtick", groups, "xticklabel", labels);
@@ -638,29 +671,47 @@
     set (gca(), "ytick", groups, "yticklabel", labels);
     hs.labels = get (gcf, "currentaxes");
   endif
+  hold off;
+
+  ## Return output arguments if desired
+  if (nargout >= 1)
+    s_o = s;
+  endif
+  if (nargout == 2)
+    hs_o = hs;
+  endif
 
 endfunction
 
-function htags = plot_tags (out_tags_x, out_tags_y, out_idx, out_IDs, sample_IDs, opt)
-  for i=1:length (out_tags_x)
+
+function htags = plot_tags (out_tags_x, out_tags_y, out_idx, out_IDs, ...
+                 sample_IDs, opt)
+
+  for i=1 : length (out_tags_x)
     if (! isempty (sample_IDs))
-      htags(i) = text(out_tags_x(i), out_tags_y(i), out_IDs{i}, opt{});
+      htags(i) = text (out_tags_x(i), out_tags_y(i), out_IDs{i}, opt{:});
     else
-      htags(i) = text(out_tags_x(i), out_tags_y(i), num2str(out_idx(i)), opt{});
+      htags(i) = text (out_tags_x(i), out_tags_y(i), num2str (out_idx(i)), ...
+                       opt{:});
     endif
   endfor
+
 endfunction
 
+
 %!demo
-%! axis ([0,3]);
-%! boxplot ({randn(10,1)*5+140, randn(13,1)*8+135});
-%! set(gca (), "xtick", [1 2], "xticklabel", {"girls", "boys"})
+%! axis ([0, 3]);
+%! boxplot ({(randn(10, 1) * 5 + 140), (randn (13, 1) * 8 + 135)});
+%! set (gca (), "xtick", [1 2], "xticklabel", {"girls", "boys"})
 %! title ("Grade 3 heights");
 
 %!demo
-%! data = [randn(10,1)*5+140; randn(25,1)*8+135; randn(20,1)*6+165];
-%! groups = [ones(10,1); ones(25,1)*2; ones(20,1)*3];
+%! data = [(randn (10, 1) * 5 + 140); (randn (25, 1) * 8 + 135); ...
+%!         (randn (20, 1) * 6 + 165)];
+%! groups = [(ones (10, 1)); (ones (25, 1) * 2); (ones (20, 1) * 3)];
 %! labels = {"Team A", "Team B", "Team C"};
-%! pos = [2,1,3];
-%! boxplot (data,groups,"Labels", labels, "Positions", pos, "OutlierTags", "on", "BoxStyle", "filled");
+%! pos = [2, 1, 3];
+%! boxplot (data, groups, "Notch", "on", "Labels", labels, "Positions", pos, "OutlierTags", ...
+%!          "on", "BoxStyle", "filled");
 %! title ("Example of Group splitting with paired vectors");
+
