A script for plotting Tukey's rootograms.

rootogram.m

% [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