classdef (InferiorClasses = {?matlab.graphics.axis.Axes, ?matlab.ui.control.UIAxes}) UnivariateDistribution < prob.ProbabilityDistribution
%UnivariateDistribution Univariate probability distribution.

%   Copyright 2012-2023 The MathWorks, Inc.

    methods(Abstract = true)
        % Method added for univariate class
        [varargout] = icdf(varargin)
    end

    methods
        % Simple implementation often overridden in derived class
        function y = random(this,varargin)
%RANDOM Random number generation.
%    R = RANDOM(P) generates a random number R from the probability
%    distribution P.
%
%    R = RANDOM(P,M,N,...) or R = RANDOM(P,[M,N,...]) creates an
%    M-by-N-by-... array R of independent random numbers from P.
%
%    Example: Generate normal random variables with mean 100 and standard
%             deviation 10.
%        p = makedist('normal', 'mu',100, 'sigma',10)
%        x = random(p,1000,1);
%        hist(x,40)
%
%    See also FITDIST, MAKEDIST.
            requireScalar(this);
            y = icdf(this,rand(varargin{:}));
        end

        % Methods suitable for univariate classes
        function y = median(this)
%MEDIAN Median.
%    M = MEDIAN(P) computes the median M for the probability distribution
%    P. M is the value such that half of the probabilty is below M and half
%    is above.
%
%    Example: Create a Weibull distribution and demonstrate that its mean
%             is larger than its median.
%        p = makedist('weibull', 'A',200, 'B',1.1)
%        wmean = mean(p)
%        wmedian = median(p)
%
%    See also IQR, ICDF.
            y = icdf(this,.5);
        end
        function y = iqr(this)
%IQR Interquartile range.
%    R = IQR(P) computes the interquartile range R for the probability
%    distribution P. R is the difference between the 75th and 25th
%    percentage points of P.
%
%    Example: Compute the interquartile range for the normal and t
%             distributions with the same location and scale parameters.
%        p = makedist('normal', 'mu',100, 'sigma',10)
%        niqr = iqr(p)
%        q = makedist('tlocationscale', 'mu',100, 'sigma',10, 'nu',5)
%        tiqr = iqr(q)
%
%    See also MEDIAN, ICDF.
            y = icdf(this,.75) - icdf(this,.25);
        end

        function hout = plot(varargin)
%PLOT Plot probability distribution.
%    PLOT(P) plots a pdf of the probability distribution P. If P was
%    created by fitting to data, the pdf is superimposed over a histogram
%    of the data.
%
%    PLOT(P,'NAME1',VALUE1,'NAME2,Value2,...) specifies optional name/value
%    arguments.
%
%      'PlotType'  - Plot type, specified as one of the following:
%                     'pdf'         Probability density function,
%                                   superimposed on a histogram if P was
%                                   created by fitting to data. This is the
%                                   default.
%                     'cdf'         Cumulative distribution function,
%                                   superimposed on an empirical cdf if P
%                                   was created by fitting to data.
%                     'probability' Probability plot. Requires that P was
%                                   created by fitting to data.
%      'Discrete'  - Indicator that specifies whether to plot the density
%                    as a discrete distribution or as a continuous
%                    distribution. When the value is true, the 'pdf' plot
%                    is a stem plot and the 'cdf' plot is a stairs plot.
%                    This option is used only if P represents a discrete
%                    distribution. Default is true.
%      'Parent'    - Handle AX for the axes to receive the plot. Default
%                    is GCA.
%
%    PLOT(AX,...) plots into the axes specified by the handle AX.
%
%    H = PLOT(...) returns handles to the plotted objects.
%
%    Example:
%        load carsmall Weight
%        p = fitdist(Weight,'lognormal')
%        plot(p)

            [varargin{:}] = convertStringsToChars(varargin{:});

            % Check for axes
            [ax,varargin] = axescheck(varargin{:});
            p = varargin{1}; % cannot get to this method without at least one element
            varargin(1) = [];

            % If SMLT not available, only pdf allowed
            args  = {'PlotType', 'Parent' 'Discrete'};
            dflts = {'pdf',      ax       true};
            [ptype,ax,isdiscrete] = statslib.internal.parseArgs(args,dflts,varargin{:});

            if ~isempty(ax) && ~isequal(ax.Type,'axes')
                error(message('stats:probdists:BadHandle'));
            end

            alltypes = {'pdf' 'histogram' 'probability' 'cdf'};
            ptype = validatestring(ptype,alltypes,'plot','PlotType');


            if ~isempty(isdiscrete)
                if islogical(isdiscrete)
                    validateattributes(isdiscrete,{'logical'},...
                        {'scalar'},...
                        'plot','Discrete',2);
                else
                    validateattributes(isdiscrete,{'numeric'},...
                        {'scalar','integer','>=',0,'<=',1},...
                        'plot','Discrete',2);
                end
            end

            if isa(p,'prob.FittableDistribution')
                haveData = ~isempty(p.InputData);
                iscont = p.Support.iscontinuous;
            else
                haveData = false;
                iscont = ~isa(p,'prob.MultinomialDistribution');
            end

            plim = icdf(p,[.001 .999]);
                
            if haveData
                inputdata = p.InputData;
                data = inputdata.data;
                lims = [min(min(data),plim(1)), max(max(data),plim(2))];
            else
                inputdata = [];
                lims = plim;
            end

            if iscont
                xgrid = linspace(lims(1),lims(2),1001)';
            elseif isequal(isdiscrete,true) || (isempty(isdiscrete) && lims(2)-lims(1)<=201)
                iscont = false; % treat as discrete
                xgrid = (floor(lims(1)):ceil(lims(2)))';
            else
                iscont = true; % treat as continuous
                xgrid = floor(linspace(lims(1),lims(2),1001))';
            end

            switch(ptype)
                case 'cdf'
                    h = plotCDF(ax,haveData,inputdata,p,xgrid,iscont);
                case 'probability'
                    h = plotProbability(ax,haveData,inputdata,p);
                case {'pdf' 'histogram'}
                    h = plotPDF(ax,haveData,inputdata,p,xgrid,iscont);
            end
            if nargout>0
                hout = h;
            end
        end

    end
end % classdef

% -----------------------
function h = plotCDF(ax,haveData,inputdata,p,xgrid,iscont)
if isempty(ax)
    ax = gca;
end
washold = ishold(ax);

if haveData
    data = inputdata.data;
    freq = inputdata.freq;
    cens = inputdata.cens;
    [f,x] = ecdf(data,'Frequency',freq,'Censoring',cens);
    h1 = stairs(ax,x,f,'LineWidth',1);
    if ~washold
        cleanupObj = onCleanup(@() hold(ax,'off'));
        hold(ax,'on');
    end
else
    h1 = [];
end
F = cdf(p,xgrid);
if iscont
    h2 = plot(ax,xgrid,F,'LineWidth',1);
else
    h2 = stairs(ax,[xgrid(1);xgrid],[0;F],'LineWidth',1);
    ax.XLim = [min(ax.XLim(1),xgrid(1)-0.5),...
               max(ax.XLim(end),xgrid(end)+0.5)];
end
h = [h2 h1];

xlabel(ax,getString(message('stats:probdists:XAxisData')));
ylabel(ax,getString(message('stats:probdists:YAxisCDF')));
end

% -----------------------
function h = plotProbability(ax,haveData,inputdata,p)
if ~haveData
    error(message('stats:probdists:ProbPlotNoData'));
elseif ~isa(p,'prob.ParametricDistribution')
    error(message('stats:probdists:ProbPlotNoParam'));
end
data = inputdata.data;
freq = inputdata.freq;
cens = inputdata.cens;
if isempty(ax)
    h = probplot(p,data,cens,freq);
else
    h = probplot(ax,p,data,cens,freq);
end

ax = ancestor(h(1),'axes');
ax.Title.String = '';
end

% -----------------------
function h = plotPDF(ax,haveData,inputdata,p,xgrid,iscont)
if isempty(ax)
    ax = gca;
end
washold = ishold(ax);
if ~washold
    newplot(ax)
end
coi = ax.ColorOrderIndex;

if haveData
    data = inputdata.data;
    freq = inputdata.freq;
    cens = inputdata.cens;
    if isempty(freq) && isempty(cens)
        if iscont
            h1 = histogram(ax,data,'Normalization','pdf');
        else
            h1 = histogram(ax,data,'Normalization','pdf','BinMethod','integers');
        end
    else
        if isempty(cens)
            N = numel(data);
        else
            N = sum(~cens);
        end
        [F,X] = ecdf(data,'Frequency',freq,'Censoring',cens);

        if iscont
            nbins = ceil(1 + log2(N)); % Sturges rule
            [heights,centers] = ecdfhist(F,X,nbins);
        else
            centers = min(X):max(X);
            [heights,centers] = ecdfhist(F,X,centers);
        end
        if isscalar(centers)
            dx = 0.5;
        else
            dx = centers(2)-centers(1);
        end
        [~,autocolor] = matlab.graphics.chart.internal.nextstyle(ax,true,false,true);
        h1 = bar(ax,centers,heights,[centers(1)-dx,centers(end)+dx],'hist');

        % Mimic histogram by setting some properties and incrementing the
        % color order index
        h1.FaceColor = autocolor;
        h1.EdgeColor = [0 0 0];
        h1.FaceAlpha = 0.6;
        ax.ColorOrderIndex = coi+1;
    end
    if ~washold
        cleanupObj = onCleanup(@() hold(ax,'off'));
        hold(ax,'on');
    end
else
    h1 = [];
end
if iscont
    xgrid = adjustXGrid(p,xgrid);
    h2 = plot(ax,xgrid,pdf(p,xgrid),'LineWidth',1);
    adjustAxisLimits(ax);
else
    h2 = stem(ax,xgrid,pdf(p,xgrid),'LineWidth',1);
    ax.XLim = [min(ax.XLim(1),xgrid(1)-0.5),...
               max(ax.XLim(end),xgrid(end)+0.5)];
end
h = [h2 h1];

xlabel(ax,getString(message('stats:probdists:XAxisData')));
if iscont
    ylabel(ax,getString(message('stats:probdists:YAxisPDF')));
else
    ylabel(ax,getString(message('stats:probdists:YAxisProb')));
end

end

% -----------------------------------
function xgrid = adjustXGrid(p,xgrid)
% See if there appear to be finite limits on the support
ptest = icdf(p,[0 1]);
if ~any(isfinite(ptest))
    return
end

% Try to capture behavior just around the finite boundaries.
if isfinite(ptest(1))
    dp = eps(ptest(1));
    xgrid = [ptest(1)+[-dp;0;dp]; xgrid];
end
if isfinite(ptest(2))
    dp = eps(ptest(2));
    xgrid = [xgrid; ptest(2)+[-dp;0;dp]];
end
xgrid = sort(xgrid);
end

% -----------------------------------
function adjustAxisLimits(ax)
% Make sure the limits are in view
ax.XLimitMethod = 'padded';
ax.YLimitMethod = 'padded';
ax.YLim(1) = 0;
end
