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