myboxplot

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

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