Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
A script for plotting Tukey's rootograms.
% [INPUT]
% ax = An object of type "matlab.graphics.axis.Axes" representing the axis on which the rootogram is plotted (optional, default=gca).
% occ = A numeric vector representing the number of occurrences.
% exd = A numeric vector representing the expected (fitted) frequencies.
% obs = A numeric vector representing the observed (real) frequencies.
% flo = A boolean that indicates whether to use floating bars (optional, default=true).
% squ = A boolean that indicates whether to use the squared root of frequencies (optional, default=true).
%
% [OUTPUT]
% h_rec = A vector containing the rectangle handles (optional).
% h_cur = The curve handle (optional).
% h_lin = The reference line handle (optional).
%
% [NOTES]
% Inspired by the following R library: https://github.com/cran/vcd
function [h_rec,h_cur,h_lin] = rootogram(varargin)
persistent ip;
if (isempty(ip))
ip = inputParser();
ip.addRequired('occ',@(x)validateattributes(x,{'numeric'},{'vector','finite','nonempty','nonnan','real','increasing','>=',0}));
ip.addRequired('exd',@(x)validateattributes(x,{'numeric'},{'vector','finite','nonempty','nonnan','real','>=',0}));
ip.addRequired('obs',@(x)validateattributes(x,{'numeric'},{'vector','finite','nonempty','nonnan','real','>=',0}));
ip.addOptional('flo',true,@(x)validateattributes(x,{'logical'},{'scalar'}));
ip.addOptional('squ',true,@(x)validateattributes(x,{'logical'},{'scalar'}));
end
[ax,args] = axescheck(varargin{:});
ip.parse(args{:});
ip_res = ip.Results;
occ = ip_res.occ;
exd = ip_res.exd;
obs = ip_res.obs;
flo = ip_res.flo;
squ = ip_res.squ;
if (isempty(ax))
figure();
ax = gca;
end
if (~isvalid(ax))
error('Invalid axis handle specified.');
end
if (occ(1) ~= 0)
error('The first occurrence must be equal to 0');
end
occ_len = length(occ);
if (length(exd) ~= occ_len)
error('The number of expected frequencies must match the number of occurrences.');
end
if (length(obs) ~= occ_len)
error('The number of observed frequencies must match the number of occurrences.');
end
switch (nargout)
case 0
rootogram_internal(ax,occ,occ_len,exd,obs,flo,squ);
case 1
h_rec = rootogram_internal(ax,occ,occ_len,exd,obs,flo,squ);
case 2
[h_rec,h_cur] = rootogram_internal(ax,occ,occ_len,exd,obs,flo,squ);
otherwise
[h_rec,h_cur,h_lin] = rootogram_internal(ax,occ,occ_len,exd,obs,flo,squ);
end
end
function [h_rec,h_cur,h_lin] = rootogram_internal(ax,occ,occ_len,exd,obs,flo,squ)
obs_nzw = obs > 0;
if (squ)
exd = sqrt(exd);
obs = sqrt(obs);
end
b = bar(ax,obs);
b_wid = get(b,'BarWidth');
b_wid_hal = b_wid / 2;
delete(b);
rec_x = occ - b_wid_hal;
if (flo)
rec_y = exd - obs;
else
rec_y = zeros(occ_len,1);
end
hold on;
if (nargout >= 1)
h_rec = cell(occ_len,1);
for i = 1:occ_len
h_rec{i} = rectangle('Position',[rec_x(i) rec_y(i) b_wid obs(i)],'FaceColor',[0.82745 0.82745 0.82745]);
end
else
for i = 1:occ_len
rectangle('Position',[rec_x(i) rec_y(i) b_wid obs(i)],'FaceColor',[0.82745 0.82745 0.82745]);
end
end
if (nargout >= 2)
h_cur = plot(ax,occ,exd,'-r','Marker','o','MarkerFaceColor','red');
else
plot(ax,occ,exd,'-r','Marker','o','MarkerFaceColor','red');
end
ref_lin = refline(ax,[0 0]);
ref_lin.Color = 'k';
if (nargout == 3)
h_lin = ref_lin;
end
hold off;
xlabel(ax,'Number of Occurrences');
if (squ)
ylabel(ax,'Squared Root of Frequency');
else
ylabel(ax,'Frequency');
end
x_lbls = repmat({''},occ_len,1);
x_lbls(obs_nzw) = cellstr(num2str(occ(obs_nzw)));
y_tcks = get(ax,'YTick');
y_tcks_neg = y_tcks < 0;
y_tcks(y_tcks_neg) = [];
y_lbls = get(ax,'YTickLabel');
y_lbls(y_tcks_neg) = [];
set(ax,'XLim',[-1 (occ(end) + 1)],'XTick',occ,'XTickLabel',x_lbls,'YLim',[min(rec_y) (max([exd; rec_y]) * 1.05)],'YTick',y_tcks,'YTickLabel',y_lbls);
end
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.