myboxplot

  sout = myboxplot (data,notched,symbol,vertical,maxwhisker)

0001 function sout = myboxplot (data,notched,symbol,vertical,maxwhisker)
0002 %  sout = myboxplot (data,notched,symbol,vertical,maxwhisker)
0003 
0004 % Copyright (C) 2012 Dynare Team
0005 %
0006 % This file is part of Dynare.
0007 %
0008 % Dynare is free software: you can redistribute it and/or modify
0009 % it under the terms of the GNU General Public License as published by
0010 % the Free Software Foundation, either version 3 of the License, or
0011 % (at your option) any later version.
0012 %
0013 % Dynare is distributed in the hope that it will be useful,
0014 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0015 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0016 % GNU General Public License for more details.
0017 %
0018 % You should have received a copy of the GNU General Public License
0019 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0020 
0021 % % % % endif
0022 if nargin < 5 | isempty(maxwhisker), maxwhisker = 1.5; end
0023 if nargin < 4 | isempty(vertical), vertical = 1; end
0024 if nargin < 3 | isempty(symbol), symbol = ['+','o']; end
0025 if nargin < 2 | isempty(notched), notched = 0; end
0026 
0027 if length(symbol)==1, symbol(2)=symbol(1); end
0028 
0029 if notched==1, notched=0.25; end
0030 a=1-notched;
0031 
0032 % ## figure out how many data sets we have
0033 if iscell(data), 
0034   nc = length(data);
0035 else
0036 %   if isvector(data), data = data(:); end
0037   nc = size(data,2);
0038 end
0039 
0040 % ## compute statistics
0041 % ## s will contain
0042 % ##    1,5    min and max
0043 % ##    2,3,4  1st, 2nd and 3rd quartile
0044 % ##    6,7    lower and upper confidence intervals for median
0045 s = zeros(7,nc);
0046 box = zeros(1,nc);
0047 whisker_x = ones(2,1)*[1:nc,1:nc];
0048 whisker_y = zeros(2,2*nc);
0049 outliers_x = [];
0050 outliers_y = [];
0051 outliers2_x = [];
0052 outliers2_y = [];
0053 
0054 for i=1:nc
0055     %   ## Get the next data set from the array or cell array
0056     if iscell(data)
0057         col = data{i}(:);
0058     else
0059         col = data(:,i);
0060     end
0061 %   ## Skip missing data
0062 % % % % % % %   col(isnan(col) | isna (col)) = [];
0063   col(isnan(col)) = [];
0064 
0065   %   ## Remember the data length
0066   nd = length(col);
0067   box(i) = nd;
0068   if (nd > 1)
0069 %     ## min,max and quartiles
0070 %     s(1:5,i) = statistics(col)(1:5);
0071 s(1,i)=min(col);
0072 s(5,i)=max(col);
0073 s(2,i)=myprctilecol(col,25);
0074 s(3,i)=myprctilecol(col,50);
0075 s(4,i)=myprctilecol(col,75);
0076 
0077 
0078 
0079 
0080 
0081 
0082 
0083 
0084 %     ## confidence interval for the median
0085     est = 1.57*(s(4,i)-s(2,i))/sqrt(nd);
0086     s(6,i) = max([s(3,i)-est, s(2,i)]);
0087     s(7,i) = min([s(3,i)+est, s(4,i)]);
0088 %     ## whiskers out to the last point within the desired inter-quartile range
0089     IQR = maxwhisker*(s(4,i)-s(2,i));
0090     whisker_y(:,i) = [min(col(col >= s(2,i)-IQR)); s(2,i)];
0091     whisker_y(:,nc+i) = [max(col(col <= s(4,i)+IQR)); s(4,i)];
0092 %     ## outliers beyond 1 and 2 inter-quartile ranges
0093     outliers = col((col < s(2,i)-IQR & col >= s(2,i)-2*IQR) | (col > s(4,i)+IQR & col <= s(4,i)+2*IQR));
0094     outliers2 = col(col < s(2,i)-2*IQR | col > s(4,i)+2*IQR);
0095     outliers_x = [outliers_x; i*ones(size(outliers))];
0096     outliers_y = [outliers_y; outliers];
0097     outliers2_x = [outliers2_x; i*ones(size(outliers2))];
0098     outliers2_y = [outliers2_y; outliers2];
0099   elseif (nd == 1)
0100 %     ## all statistics collapse to the value of the point
0101     s(:,i) = col;
0102 %     ## single point data sets are plotted as outliers.
0103     outliers_x = [outliers_x; i];
0104     outliers_y = [outliers_y; col];
0105   else
0106 %     ## no statistics if no points
0107     s(:,i) = NaN;
0108   end
0109 end
0110 % % % % if isempty(outliers2_y)
0111 % % % %     outliers2_y=
0112 % ## Note which boxes don't have enough stats
0113 chop = find(box <= 1);
0114     
0115 % ## Draw a box around the quartiles, with width proportional to the number of
0116 % ## items in the box. Draw notches if desired.
0117 box = box*0.23/max(box);
0118 quartile_x = ones(11,1)*[1:nc] + [-a;-1;-1;1;1;a;1;1;-1;-1;-a]*box;
0119 quartile_y = s([3,7,4,4,7,3,6,2,2,6,3],:);
0120 
0121 % ## Draw a line through the median
0122 median_x = ones(2,1)*[1:nc] + [-a;+a]*box;
0123 % median_x=median(col);
0124 median_y = s([3,3],:);
0125 
0126 % ## Chop all boxes which don't have enough stats
0127 quartile_x(:,chop) = [];
0128 quartile_y(:,chop) = [];
0129 whisker_x(:,[chop,chop+nc]) = [];
0130 whisker_y(:,[chop,chop+nc]) = [];
0131 median_x(:,chop) = [];
0132 median_y(:,chop) = [];
0133 % % % %
0134 % ## Add caps to the remaining whiskers
0135 cap_x = whisker_x;
0136 cap_x(1,:) =cap_x(1,:)- 0.05;
0137 cap_x(2,:) =cap_x(2,:)+ 0.05;
0138 cap_y = whisker_y([1,1],:);
0139 
0140 % #quartile_x,quartile_y
0141 % #whisker_x,whisker_y
0142 % #median_x,median_y
0143 % #cap_x,cap_y
0144 %
0145 % ## Do the plot
0146 
0147 mm=min(min(data));
0148 MM=max(max(data));
0149 
0150 if vertical
0151     plot (quartile_x, quartile_y, 'b',  ...
0152         whisker_x, whisker_y, 'b--',   ...
0153         cap_x, cap_y, 'k',   ...
0154         median_x, median_y, 'r',  ...
0155         outliers_x, outliers_y, [symbol(1),'r'],   ...
0156         outliers2_x, outliers2_y, [symbol(2),'r']);
0157         set(gca,'XTick',1:nc);
0158         set(gca, 'XLim', [0.5, nc+0.5]);
0159         set(gca, 'YLim', [mm-(MM-mm)*0.05-eps, MM+(MM-mm)*0.05+eps]);
0160 
0161 else
0162 % % % % %     plot (quartile_y, quartile_x, "b;;",
0163 % % % % %     whisker_y, whisker_x, "b;;",
0164 % % % % %     cap_y, cap_x, "b;;",
0165 % % % % %     median_y, median_x, "r;;",
0166 % % % % %     outliers_y, outliers_x, [symbol(1),"r;;"],
0167 % % % % %     outliers2_y, outliers2_x, [symbol(2),"r;;"]);
0168 end
0169 
0170 if nargout,
0171   sout=s;
0172 end
0173 % % % endfunction