From 9efb7847632b8f805c805117622273a74fe27934 Mon Sep 17 00:00:00 2001 From: Johannes Pfeifer Date: Sat, 2 Dec 2023 12:52:01 +0100 Subject: [PATCH] GSA: remove some unused functions and move other to inline ones --- license.txt | 6 +- matlab/gsa/cumplot.m | 13 +-- matlab/gsa/gsa_plotmatrix.m | 99 -------------------- matlab/gsa/map_ident_.m | 10 ++ matlab/gsa/myboxplot.m | 113 ++++++++++++----------- matlab/gsa/myprctilecol.m | 43 --------- matlab/gsa/priorcdf.m | 5 +- matlab/gsa/read_data.m | 44 --------- matlab/gsa/stand_.m | 24 +++-- matlab/gsa/tcrit.m | 2 +- matlab/gsa/teff.m | 6 -- matlab/gsa/trank.m | 34 ------- matlab/myboxplot.m | 177 ------------------------------------ 13 files changed, 92 insertions(+), 484 deletions(-) delete mode 100644 matlab/gsa/gsa_plotmatrix.m delete mode 100644 matlab/gsa/myprctilecol.m delete mode 100644 matlab/gsa/read_data.m delete mode 100644 matlab/gsa/trank.m delete mode 100644 matlab/myboxplot.m diff --git a/license.txt b/license.txt index 71083f2ab..aa420ab47 100644 --- a/license.txt +++ b/license.txt @@ -171,16 +171,13 @@ Comment: Written by Jessica Cariboni and Francesca Campolongo Files: matlab/gsa/cumplot.m matlab/gsa/filt_mc_.m - matlab/gsa/gsa_plotmatrix.m matlab/gsa/gsa_skewness.m matlab/gsa/log_trans_.m matlab/gsa/map_calibration.m matlab/gsa/map_ident_.m matlab/gsa/mcf_analysis.m matlab/gsa/myboxplot.m - matlab/gsa/myprctilecol.m matlab/gsa/prior_draw_gsa.m - matlab/gsa/read_data.m matlab/gsa/redform_map.m matlab/gsa/redform_screen.m matlab/gsa/scatter_mcf.m @@ -191,9 +188,8 @@ Files: matlab/gsa/cumplot.m matlab/gsa/stand_.m matlab/gsa/tcrit.m matlab/gsa/teff.m - matlab/gsa/trank.m Copyright: 2011-2018 European Commission - 2011-2018 Dynare Team + 2011-2023 Dynare Team License: GPL-3+ Files: matlab/gsa/pick.m diff --git a/matlab/gsa/cumplot.m b/matlab/gsa/cumplot.m index 9de34e86c..c58567a4b 100644 --- a/matlab/gsa/cumplot.m +++ b/matlab/gsa/cumplot.m @@ -1,5 +1,10 @@ function h = cumplot(x) %function h =cumplot(x) +% Inputs: +% - x [double] data series +% +% Outputs: +% - h [handle] figure handle % Written by Marco Ratto % Joint Research Centre, The European Commission, @@ -26,9 +31,5 @@ function h = cumplot(x) n=length(x); x=[-inf; sort(x); Inf]; y=[0:n n]./n; -h0 = stairs(x,y); -grid on, - -if nargout - h=h0; -end +h = stairs(x,y); +grid on \ No newline at end of file diff --git a/matlab/gsa/gsa_plotmatrix.m b/matlab/gsa/gsa_plotmatrix.m deleted file mode 100644 index 0b5caa6bf..000000000 --- a/matlab/gsa/gsa_plotmatrix.m +++ /dev/null @@ -1,99 +0,0 @@ -function gsa_plotmatrix(type,varargin) -% function gsa_plotmatrix(type,varargin) -% extended version of the standard MATLAB plotmatrix -% -% Written by Marco Ratto -% Joint Research Centre, The European Commission, -% marco.ratto@ec.europa.eu - -% Copyright © 2011-2012 European Commission -% Copyright © 2011-2017 Dynare Team -% -% This file is part of Dynare. -% -% Dynare is free software: you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation, either version 3 of the License, or -% (at your option) any later version. -% -% Dynare is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with Dynare. If not, see . - -global bayestopt_ options_ M_ - -RootDirectoryName = CheckPath('gsa',M_.dname); - -if options_.opt_gsa.pprior - load([ RootDirectoryName filesep M_.fname '_prior.mat'],'lpmat0','lpmat','istable','iunstable','iindeterm','iwrong') -else - load([ RootDirectoryName filesep M_.fname '_mc.mat'],'lpmat0','lpmat','istable','iunstable','iindeterm','iwrong') - eval(['load ' options_.mode_file ' xparam1;']'); -end - -iexplosive = iunstable(~ismember(iunstable,[iindeterm;iwrong])); - -switch type - case 'all' - x=[lpmat0 lpmat]; - NumberOfDraws=size(x,1); - B=NumberOfDraws; - case 'stable' - x=[lpmat0(istable,:) lpmat(istable,:)]; - NumberOfDraws=size(x,1); - B=NumberOfDraws; - case 'nosolution' - x=[lpmat0(iunstable,:) lpmat(iunstable,:)]; - NumberOfDraws=size(x,1); - B=NumberOfDraws; - case 'unstable' - x=[lpmat0(iexplosive,:) lpmat(iexplosive,:)]; - NumberOfDraws=size(x,1); - B=NumberOfDraws; - case 'indeterm' - x=[lpmat0(iindeterm,:) lpmat(iindeterm,:)]; - NumberOfDraws=size(x,1); - B=NumberOfDraws; - case 'wrong' - x=[lpmat0(iwrong,:) lpmat(iwrong,:)]; - NumberOfDraws=size(x,1); - B=NumberOfDraws; - -end - -if isempty(x) - disp('Empty parameter set!') - return -end - -for j=1:length(varargin) - jcol(j)=strmatch(varargin{j},bayestopt_.name,'exact'); -end - -[H,AX,BigA,P,PAx]=plotmatrix(x(:,jcol)); - -for j=1:length(varargin) - % axes(AX(1,j)), title(varargin{j}) - % axes(AX(j,1)), ylabel(varargin{j}) - % set(AX(1,j),'title',varargin{j}), - set(get(AX(j,1),'ylabel'),'string',varargin{j}) - set(get(AX(end,j),'xlabel'),'string',varargin{j}) -end - -if options_.opt_gsa.pprior==0 - xparam1=xparam1(jcol); - for j=1:length(varargin) - for i=1:j-1 - axes(AX(j,i)) - hold on, plot(xparam1(i),xparam1(j),'*r') - end - for i=j+1:length(varargin) - axes(AX(j,i)) - hold on, plot(xparam1(i),xparam1(j),'*r') - end - end -end diff --git a/matlab/gsa/map_ident_.m b/matlab/gsa/map_ident_.m index ec84adef2..2b7194fa2 100644 --- a/matlab/gsa/map_ident_.m +++ b/matlab/gsa/map_ident_.m @@ -365,3 +365,13 @@ if options_.TeX && any(strcmp('eps',cellstr(options_.graph_format))) fprintf(fidTeX,'%% End Of TeX file. \n'); fclose(fidTeX); end + + +function yr = trank(y) +% yr is the rank transformation of y +yr=NaN(size(y)); +[nr, nc] = size(y); +for j=1:nc + [~, is]=sort(y(:,j)); + yr(is,j)=[1:nr]'./nr; +end diff --git a/matlab/gsa/myboxplot.m b/matlab/gsa/myboxplot.m index 7864107e2..4d6cf60d1 100644 --- a/matlab/gsa/myboxplot.m +++ b/matlab/gsa/myboxplot.m @@ -1,12 +1,8 @@ function sout = myboxplot (data,notched,symbol,vertical,maxwhisker) -% sout = myboxplot (data,notched,symbol,vertical,maxwhisker) +% sout = myboxplot (data,notched,symbol,vertical,maxwhisker) +% Creates a box plot -% Written by Marco Ratto -% Joint Research Centre, The European Commission, -% marco.ratto@ec.europa.eu - -% Copyright © 2012 European Commission -% Copyright © 2012-2017 Dynare Team +% Copyright © 2010-2023 Dynare Team % % This file is part of Dynare. % @@ -23,18 +19,17 @@ function sout = myboxplot (data,notched,symbol,vertical,maxwhisker) % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . -% % % % endif -if nargin < 5 | isempty(maxwhisker), maxwhisker = 1.5; end -if nargin < 4 | isempty(vertical), vertical = 1; end -if nargin < 3 | isempty(symbol), symbol = ['+','o']; end -if nargin < 2 | isempty(notched), notched = 0; end +if nargin < 5 || isempty(maxwhisker), maxwhisker = 1.5; end +if nargin < 4 || isempty(vertical), vertical = 1; end +if nargin < 3 || isempty(symbol), symbol = ['+','o']; end +if nargin < 2 || isempty(notched), notched = 0; end if length(symbol)==1, symbol(2)=symbol(1); end if notched==1, notched=0.25; end a=1-notched; -% ## figure out how many data sets we have +% % figure out how many data sets we have if iscell(data) nc = length(data); else @@ -42,11 +37,11 @@ else nc = size(data,2); end -% ## 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 +% 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); whisker_x = ones(2,1)*[1:nc,1:nc]; @@ -57,44 +52,36 @@ outliers2_x = []; outliers2_y = []; for i=1:nc - % ## Get the next data set from the array or cell array + % Get the next data set from the array or cell array if iscell(data) col = data{i}(:); else col = data(:,i); end - % ## Skip missing data + % Skip missing data % % % % % % % col(isnan(col) | isna (col)) = []; col(isnan(col)) = []; - % ## Remember the data length + % Remember the data length nd = length(col); box(i) = nd; if (nd > 1) - % ## min,max and quartiles - % s(1:5,i) = statistics(col)(1:5); + % min,max and quartiles s(1,i)=min(col); s(5,i)=max(col); s(2,i)=myprctilecol(col,25); s(3,i)=myprctilecol(col,50); s(4,i)=myprctilecol(col,75); - - - - - - - - % ## confidence interval for the median + % confidence interval for the median est = 1.57*(s(4,i)-s(2,i))/sqrt(nd); s(6,i) = max([s(3,i)-est, s(2,i)]); s(7,i) = min([s(3,i)+est, s(4,i)]); - % ## whiskers out to the last point within the desired inter-quartile range + % whiskers out to the last point within the desired inter-quartile range IQR = maxwhisker*(s(4,i)-s(2,i)); whisker_y(:,i) = [min(col(col >= s(2,i)-IQR)); s(2,i)]; whisker_y(:,nc+i) = [max(col(col <= s(4,i)+IQR)); s(4,i)]; - % ## outliers beyond 1 and 2 inter-quartile ranges + % outliers beyond 1 and 2 inter-quartile ranges outliers = col((col < s(2,i)-IQR & col >= s(2,i)-2*IQR) | (col > s(4,i)+IQR & col <= s(4,i)+2*IQR)); outliers2 = col(col < s(2,i)-2*IQR | col > s(4,i)+2*IQR); outliers_x = [outliers_x; i*ones(size(outliers))]; @@ -102,41 +89,37 @@ for i=1:nc outliers2_x = [outliers2_x; i*ones(size(outliers2))]; outliers2_y = [outliers2_y; outliers2]; elseif (nd == 1) - % ## all statistics collapse to the value of the point + % all statistics collapse to the value of the point s(:,i) = col; - % ## single point data sets are plotted as outliers. + % single point data sets are plotted as outliers. outliers_x = [outliers_x; i]; outliers_y = [outliers_y; col]; else - % ## no statistics if no points + % no statistics if no points s(:,i) = NaN; end end -% % % % if isempty(outliers2_y) -% % % % outliers2_y= -% ## Note which boxes don't have enough stats +% Note which boxes don't have enough stats chop = find(box <= 1); -% ## Draw a box around the quartiles, with width proportional to the number of -% ## items in the box. Draw notches if desired. +% Draw a box around the quartiles, with width proportional to the number of +% items in the box. Draw notches if desired. box = box*0.23/max(box); quartile_x = ones(11,1)*[1:nc] + [-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 +% Draw a line through the median median_x = ones(2,1)*[1:nc] + [-a;+a]*box; -% median_x=median(col); median_y = s([3,3],:); -% ## Chop all boxes which don't have enough stats +% 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]) = []; median_x(:,chop) = []; median_y(:,chop) = []; -% % % % -% ## Add caps to the remaining whiskers +% Add caps to the remaining whiskers cap_x = whisker_x; cap_x(1,:) =cap_x(1,:)- 0.05; cap_x(2,:) =cap_x(2,:)+ 0.05; @@ -146,11 +129,14 @@ cap_y = whisker_y([1,1],:); % #whisker_x,whisker_y % #median_x,median_y % #cap_x,cap_y -% -% ## Do the plot +% Do the plot mm=min(min(data)); MM=max(max(data)); +if isnan(mm) + mm=0; + MM=0; +end if vertical plot (quartile_x, quartile_y, 'b', ... @@ -162,17 +148,30 @@ if vertical set(gca,'XTick',1:nc); set(gca, 'XLim', [0.5, nc+0.5]); set(gca, 'YLim', [mm-(MM-mm)*0.05-eps, MM+(MM-mm)*0.05+eps]); - -else - % % % % % 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;;"]); end if nargout sout=s; end -% % % endfunction \ No newline at end of file + +function y = myprctilecol(x,p) + +xx = sort(x); +[m,n] = size(x); + +if m==1 | n==1 + m = max(m,n); + if m == 1 + y = x*ones(length(p),1); + return + end + n = 1; + q = 100*(0.5:m - 0.5)./m; + xx = [min(x); xx(:); max(x)]; +else + q = 100*(0.5:m - 0.5)./m; + xx = [min(x); xx; max(x)]; +end + +q = [0 q 100]; +y = interp1(q,xx,p); \ No newline at end of file diff --git a/matlab/gsa/myprctilecol.m b/matlab/gsa/myprctilecol.m deleted file mode 100644 index 1baafb1e0..000000000 --- a/matlab/gsa/myprctilecol.m +++ /dev/null @@ -1,43 +0,0 @@ -function y = myprctilecol(x,p) - -% Written by Marco Ratto -% Joint Research Centre, The European Commission, -% marco.ratto@ec.europa.eu - -% Copyright © 2012 European Commission -% Copyright © 2012-2017 Dynare Team -% -% This file is part of Dynare. -% -% Dynare is free software: you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation, either version 3 of the License, or -% (at your option) any later version. -% -% Dynare is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with Dynare. If not, see . - -xx = sort(x); -[m,n] = size(x); - -if m==1 | n==1 - m = max(m,n); - if m == 1 - y = x*ones(length(p),1); - return - end - n = 1; - q = 100*(0.5:m - 0.5)./m; - xx = [min(x); xx(:); max(x)]; -else - q = 100*(0.5:m - 0.5)./m; - xx = [min(x); xx; max(x)]; -end - -q = [0 q 100]; -y = interp1(q,xx,p); \ No newline at end of file diff --git a/matlab/gsa/priorcdf.m b/matlab/gsa/priorcdf.m index 1be3e1f14..ac0ffc1a6 100644 --- a/matlab/gsa/priorcdf.m +++ b/matlab/gsa/priorcdf.m @@ -1,5 +1,5 @@ function xcum = priorcdf(para, pshape, p6, p7, p3, p4) - +% xcum = priorcdf(para, pshape, p6, p7, p3, p4) % This procedure transforms x vectors into cumulative values % pshape: 0 is point mass, both para and p2 are ignored % 1 is BETA(mean,stdd) @@ -11,7 +11,7 @@ function xcum = priorcdf(para, pshape, p6, p7, p3, p4) % 8 is WEIBULL(s, k) % Adapted by M. Ratto from MJ priordens.m -% Copyright © 2012-2015 Dynare Team +% Copyright © 2012-2023 Dynare Team % % This file is part of Dynare. % @@ -28,6 +28,7 @@ function xcum = priorcdf(para, pshape, p6, p7, p3, p4) % You should have received a copy of the GNU General Public License % along with Dynare. If not, see . +xcum=NaN(size(para)); for i=1:length(pshape) switch pshape(i) case 1 % (generalized) BETA Prior diff --git a/matlab/gsa/read_data.m b/matlab/gsa/read_data.m deleted file mode 100644 index 5af0f5204..000000000 --- a/matlab/gsa/read_data.m +++ /dev/null @@ -1,44 +0,0 @@ -function [gend, data] = read_data() -% Written by Marco Ratto -% Joint Research Centre, The European Commission, -% marco.ratto@ec.europa.eu - -% Copyright © 2012-2015 European Commission -% Copyright © 2012-2017 Dynare Team -% -% This file is part of Dynare. -% -% Dynare is free software: you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation, either version 3 of the License, or -% (at your option) any later version. -% -% Dynare is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with Dynare. If not, see . - -global options_ - -rawdata = read_variables(options_.datafile,options_.varobs,[],options_.xls_sheet,options_.xls_range); - -options_ = set_default_option(options_,'nobs',size(rawdata,1)-options_.first_obs+1); -gend = options_.nobs; - -rawdata = rawdata(options_.first_obs:options_.first_obs+gend-1,:); -if options_.loglinear == 1 & ~options_.logdata - rawdata = log(rawdata); -end -if options_.prefilter == 1 - data = transpose(rawdata-ones(gend,1)* mean(rawdata,1)); -else - data = transpose(rawdata); -end - -if ~isreal(rawdata) - error(['There are complex values in the data. Probably a wrong' ... - ' transformation']) -end diff --git a/matlab/gsa/stand_.m b/matlab/gsa/stand_.m index 176af7f5c..c95ccf143 100644 --- a/matlab/gsa/stand_.m +++ b/matlab/gsa/stand_.m @@ -1,20 +1,22 @@ function [y, meany, stdy] = stand_(x) -% STAND_ Standardise a matrix by columns +% [y, meany, stdy] = stand_(x) +% Standardise a matrix by columns % % [x,my,sy]=stand(y) % -% y: Time series (column matrix) +% Inputs: +% - x: Time series (column matrix) % -% x: standardised equivalent of y -% my: Vector of mean values for each column of y -% sy: Vector of standard deviations for each column of y +% - y: standardised equivalent of x +% - meany: Vector of mean values for each column of x +% - stdy: Vector of standard deviations for each column of x % % Written by Marco Ratto % Joint Research Centre, The European Commission, % marco.ratto@ec.europa.eu % Copyright © 2012 European Commission -% Copyright © 2012-2017 Dynare Team% +% Copyright © 2012-2023 Dynare Team % This file is part of Dynare. % % Dynare is free software: you can redistribute it and/or modify @@ -34,9 +36,11 @@ if nargin==0 return end +meany=NaN(size(x,2),1); +stdy=NaN(size(x,2),1); +y=NaN(size(x)); for j=1:size(x,2) - meany(j)=mean(x(find(~isnan(x(:,j))),j)); - stdy(j)=std(x(find(~isnan(x(:,j))),j)); + meany(j)=mean(x(~isnan(x(:,j)),j)); + stdy(j)=std(x(~isnan(x(:,j)),j)); y(:,j)=(x(:,j)-meany(j))./stdy(j); -end -% end of m-file \ No newline at end of file +end \ No newline at end of file diff --git a/matlab/gsa/tcrit.m b/matlab/gsa/tcrit.m index c6b17f4bc..120a6107b 100644 --- a/matlab/gsa/tcrit.m +++ b/matlab/gsa/tcrit.m @@ -150,4 +150,4 @@ if n<=100 t_crit=t_crit(n,ncol); else t_crit=t_crit(end,ncol); -end \ No newline at end of file +end diff --git a/matlab/gsa/teff.m b/matlab/gsa/teff.m index 50c4917ee..b815237bd 100644 --- a/matlab/gsa/teff.m +++ b/matlab/gsa/teff.m @@ -37,15 +37,12 @@ if ndim==3 [ir, ic]=(find( (tmax-tmin)>1.e-8)); j0 = length(ir); yt=zeros(Nsam, j0); - for j=1:j0 y0=squeeze(T(ir(j),ic(j),:)); - %y1=ones(size(lpmat,1),1)*NaN; y1=ones(Nsam,1)*NaN; y1(istable,1)=y0; yt(:,j)=y1; end - else tmax=max(T,[],2); tmin=min(T,[],2); @@ -53,7 +50,4 @@ else j0 = length(ir); yt=NaN(Nsam, j0); yt(istable,:)=T(ir,:)'; - - end -%clear y0 y1; diff --git a/matlab/gsa/trank.m b/matlab/gsa/trank.m deleted file mode 100644 index 41a200a94..000000000 --- a/matlab/gsa/trank.m +++ /dev/null @@ -1,34 +0,0 @@ -function yr = trank(y) -% yr = trank(y); -% yr is the rank transformation of y -% -% Written by Marco Ratto -% Joint Research Centre, The European Commission, -% marco.ratto@ec.europa.eu -% -% Reference: -% M. Ratto, Global Sensitivity Analysis for Macroeconomic models, MIMEO, 2006. - -% Copyright © 2012 European Commission -% Copyright © 2012-2017 Dynare Team -% -% This file is part of Dynare. -% -% Dynare is free software: you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation, either version 3 of the License, or -% (at your option) any later version. -% -% Dynare is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with Dynare. If not, see . - -[nr, nc] = size(y); -for j=1:nc - [dum, is]=sort(y(:,j)); - yr(is,j)=[1:nr]'./nr; -end diff --git a/matlab/myboxplot.m b/matlab/myboxplot.m deleted file mode 100644 index 73eb78987..000000000 --- a/matlab/myboxplot.m +++ /dev/null @@ -1,177 +0,0 @@ - -function sout = myboxplot (data,notched,symbol,vertical,maxwhisker) - -% sout = myboxplot (data,notched,symbol,vertical,maxwhisker) - -% -% Copyright © 2010-2017 Dynare Team -% -% This file is part of Dynare. -% -% Dynare is free software: you can redistribute it and/or modify -% it under the terms of the GNU General Public License as published by -% the Free Software Foundation, either version 3 of the License, or -% (at your option) any later version. -% -% Dynare is distributed in the hope that it will be useful, -% but WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -% GNU General Public License for more details. -% -% You should have received a copy of the GNU General Public License -% along with Dynare. If not, see . - -% % % % endif -if nargin < 5 || isempty(maxwhisker), maxwhisker = 1.5; end -if nargin < 4 || isempty(vertical), vertical = 1; end -if nargin < 3 || isempty(symbol), symbol = ['+','o']; end -if nargin < 2 || isempty(notched), notched = 0; end - -if length(symbol)==1, symbol(2)=symbol(1); end - -if notched==1, notched=0.25; end -a=1-notched; - -% ## figure out how many data sets we have -if iscell(data) - nc = length(data); -else - % if isvector(data), data = data(:); end - nc = size(data,2); -end - -% ## 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); -whisker_x = ones(2,1)*[1:nc,1:nc]; -whisker_y = zeros(2,2*nc); -outliers_x = []; -outliers_y = []; -outliers2_x = []; -outliers2_y = []; - -for i=1:nc - % ## Get the next data set from the array or cell array - if iscell(data) - col = data{i}(:); - else - col = data(:,i); - end - % ## Skip missing data - % % % % % % % col(isnan(col) | isna (col)) = []; - col(isnan(col)) = []; - - % ## Remember the data length - nd = length(col); - box(i) = nd; - if (nd > 1) - % ## min,max and quartiles - % s(1:5,i) = statistics(col)(1:5); - s(1,i)=min(col); - s(5,i)=max(col); - s(2,i)=myprctilecol(col,25); - s(3,i)=myprctilecol(col,50); - s(4,i)=myprctilecol(col,75); - - - - - - - - - % ## confidence interval for the median - est = 1.57*(s(4,i)-s(2,i))/sqrt(nd); - s(6,i) = max([s(3,i)-est, s(2,i)]); - s(7,i) = min([s(3,i)+est, s(4,i)]); - % ## whiskers out to the last point within the desired inter-quartile range - IQR = maxwhisker*(s(4,i)-s(2,i)); - whisker_y(:,i) = [min(col(col >= s(2,i)-IQR)); s(2,i)]; - whisker_y(:,nc+i) = [max(col(col <= s(4,i)+IQR)); s(4,i)]; - % ## outliers beyond 1 and 2 inter-quartile ranges - outliers = col((col < s(2,i)-IQR & col >= s(2,i)-2*IQR) | (col > s(4,i)+IQR & col <= s(4,i)+2*IQR)); - outliers2 = col(col < s(2,i)-2*IQR | col > s(4,i)+2*IQR); - outliers_x = [outliers_x; i*ones(size(outliers))]; - outliers_y = [outliers_y; outliers]; - outliers2_x = [outliers2_x; i*ones(size(outliers2))]; - outliers2_y = [outliers2_y; outliers2]; - elseif (nd == 1) - % ## all statistics collapse to the value of the point - s(:,i) = col; - % ## single point data sets are plotted as outliers. - outliers_x = [outliers_x; i]; - outliers_y = [outliers_y; col]; - else - % ## no statistics if no points - s(:,i) = NaN; - end -end -% % % % if isempty(outliers2_y) -% % % % outliers2_y= -% ## Note which boxes don't have enough stats -chop = find(box <= 1); - -% ## Draw a box around the quartiles, with width proportional to the number of -% ## items in the box. Draw notches if desired. -box = box*0.23/max(box); -quartile_x = ones(11,1)*[1:nc] + [-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)*[1:nc] + [-a;+a]*box; -% median_x=median(col); -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]) = []; -median_x(:,chop) = []; -median_y(:,chop) = []; -% % % % -% ## Add caps to the remaining whiskers -cap_x = whisker_x; -cap_x(1,:) =cap_x(1,:)- 0.05; -cap_x(2,:) =cap_x(2,:)+ 0.05; -cap_y = whisker_y([1,1],:); - -% #quartile_x,quartile_y -% #whisker_x,whisker_y -% #median_x,median_y -% #cap_x,cap_y -% -% ## Do the plot - -mm=min(min(data)); -MM=max(max(data)); -if isnan(mm), mm=0; MM=0; end - -if vertical - plot (quartile_x, quartile_y, 'b', ... - whisker_x, whisker_y, 'b--', ... - cap_x, cap_y, 'k', ... - median_x, median_y, 'r', ... - outliers_x, outliers_y, [symbol(1),'r'], ... - outliers2_x, outliers2_y, [symbol(2),'r']); - set(gca,'XTick',1:nc); - set(gca, 'XLim', [0.5, nc+0.5]); - set(gca, 'YLim', [mm-(MM-mm)*0.05-eps, MM+(MM-mm)*0.05+eps]); - -else - % % % % % 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;;"]); -end - -if nargout - sout=s; -end -% % % endfunction \ No newline at end of file