classdef (InferiorClasses = {?matlab.graphics.axis.Axes, ?matlab.ui.control.UIAxes, ?matlab.ui.Figure}) ...
        MultinomialRegression < classreg.regr.ParametricRegression
    %MULTINOMIALREGRESSION Multinomial regression for multi-class classification.
    %   MultinomialRegression is an object storing a fitted multinomial
    %   regression to be used for binary or multi-class classification
    %   problems.
    %
    %   Create a MultinomialRegression object using the fitmnr method, and
    %   providing data in one of the following ways:
    %       - Provide a double/single matrix, X, of predictor data, and
    %         provide a vector of classes Y. Y can also be an N x K matrix
    %         of counts, where there are N observations, K classes, and
    %         Y(i,j) is the number of times observation i has response j.
    %       - Provide a table of predictor data, X, and Y as specified
    %         above.
    %       - Provide a table of data, TBL, and the name of a column of TBL
    %         to treat as the response, or a formula in Wilkinson notation
    %         using the names in TBL as the terms in the formula.
    %
    %   MULTINOMIALREGRESSION Properties
    %   ClassNames                - Names of classes used to fit the model
    %   CoefficientCovariance     - Covariance matrix of the estimated coefficients
    %   CoefficientNames          - Names of the estimated coefficients
    %   Coefficients              - Table of coefficients and relevant statistics
    %   Deviance                  - Deviance of the fitted model
    %   DFE                       - Degrees of freedom of error of the model
    %   Dispersion                - Estimated or theoretical dispersion value used by the model
    %   EstimateDispersion        - Indicator of if the dispersion value is estimated in fitting
    %   Fitted                    - Predicted response values from the training data
    %   Formula                   - Formula representing the model specification
    %   IncludeClassInteractions  - Indicator of if the model was fit with interactions between
    %                               the coefficients and classes
    %   Link                      - Link function for the model
    %   Loglikelihood             - Log-likelihood of the fitted model
    %   ModelCriterion            - Struct of a few model fit criteria
    %   ModelType                 - One of "nominal", "ordinal", or "hierarchical", indicating
    %                               what kind of model the object represents
    %   NumCoefficients           - Number of coefficients in the model
    %   NumObservations           - Number of observations used to fit the model
    %   NumPredictors             - Number of predictors used to fit the model
    %   NumVariables              - Number of variables used in creating the model
    %   ObservationInfo           - Table of information about the observations, including
    %                               missing observations and observation weights
    %   ObservationNames          - Names of the observations used to fit the model
    %   PredictorNames            - Names of the predictors used in the model
    %   Residuals                 - Table of a few fitted model residuals
    %   ResponseName              - Name of the response variable for the model
    %   Rsquared                  - Struct of two kinds of pseudo R^2 (coefficient of determination)
    %                               values of the fitted modelE
    %   VariableInfo              - Table of information about the data used to fit the
    %                               model, including ranges of data
    %   VariableNames             - Names of the variables used to fit the model
    %   Variables                 - Table of the variables used to fit the model
    %
    %   MULTINOMIALREGRESSION Methods
    %   coefCI                - Get confidence interval estimates for model coefficients
    %   coefTest              - Run linear hypothesis testing on the model coefficients
    %   feval                 - Predict responses given predictor data in separate variables
    %   partialDependence     - Compute the partial dependence between the
    %                           predictors and response
    %   plotPartialDependence - Plot the partial dependence between the
    %                           predictors and response
    %   predict               - Predict responses, probabilities, and confidence bounds
    %                           for a set of predictor data
    %   random                - Generate random samples of responses given a set of
    %                           predictor data
    %   testDeviance          - Test the goodness-of-fit of the model using the deviance
    %   plotSlice             - Generate different styles of plots through the fitted
    %                           regression surface
    %   plotResiduals         - Generate different styles of plots of the model residuals
    %
    %   Example: Fit a multinomial regression to the Fisher iris dataset
    %   load fisheriris
    %   mdl = fitmnr(meas, species)

    %   Copyright 2022-2023 The MathWorks, Inc.

    properties (GetAccess=public, SetAccess=protected)
        %FITTED - A vector of fitted (predicted) response values for the
        %   training data used to fit the MultinomialRegression.
        Fitted

        %MODELTYPE - One of 'nominal', 'ordinal', or 'hierarchical',
        %   indicating which kind of model the MultinomialRegression
        %   represents.
        ModelType

        %INCLUDECLASSINTERACTIONS - A logical scalar indicating if the
        %   MultinomialRegression object was fit with interactions between
        %   the coefficients and the classes. If false, the model has a
        %   common set of coefficients for the predictor variables across
        %   all multinomial categories.  This is often described as "parallel
        %   regression". If true, the model has different coefficients
        %   across categories. In all cases, the model has different
        %   intercepts across categories.
        IncludeClassInteractions

        %DEVIANCE - A double/single scalar indicating the deviance of the
        %   fitted MultinomialRegression. It is twice the difference between the
        %   maximum achievable log likelihood and that attained by the fitted model.
        Deviance

        %RESIDUALS - A table of the residuals of the fitted MultinomialRegression.
        %   It has one row per observation, and columns for the following
        %   residual types:
        %       'Raw'      - Observed minus fitted values
        %       'Pearson'  - Raw residuals divided by the estimated standard
        %                    deviation of the response.
        %       'Deviance' - Residuals based on the contribution of each
        %                    observation to the deviance.
        Residuals
    end

    properties (Dependent, GetAccess=public, SetAccess=private)
        %ESTIMATEDISPERSION - A logical scalar indicating if the fitted
        %   MultinomialRegression estimated a dispersion parameter, or used
        %   the theoretical dispersion value of 1.
        EstimateDispersion

        %DISPERSION - A double/single scalar of the dispersion of the
        %   MultinomialRegression. If EstimateDispersion is true, this is
        %   the estimated dispersion parameter. If false, this is the
        %   theoretical value of 1.
        Dispersion

        %LINK - One of 'logit', 'probit', 'comploglog', 'loglog',
        %   indicating what kind of link function was used for the
        %   MultinomialRegression.
        Link

        %CLASSNAMES - A vector of the same type as Fitted. It contains the
        %   unique set of classes used to fit the MultinomialRegression.
        ClassNames
    end

    properties (GetAccess=private, SetAccess=private)
        % PrivFitStats - A struct that stores the stats output from mnrfit.
        % Used to populate some properties, and to get outputs for some
        % functions
        PrivFitStats

        % CoefSize - A 2 element vector that contains the original size of
        % the fitted coefficients. The object treats the coefficients as a
        % single vector, but mnrval expects a matrix
        CoefSize

        % PrivClassNames - A classreg.learning.internal.ClassLabel object
        % storing the class names. Used to populate model.ClassNames
        PrivClassNames
    end

    methods
        %%% GET METHODS %%%
        function ed = get.EstimateDispersion(model)
            ed = model.PrivFitStats.estdisp;
        end

        function d = get.Dispersion(model)
            % mnrfit treats the dispersion as the standard deviation. Most
            % other texts on GLMs treat dispersion as the variance. Treat
            % it as the variance here to be consistent with GeneralizedLinearModel
            d = model.PrivFitStats.s^2;
        end

        function link = get.Link(model)
            link = model.Formula.Link;
        end

        function cl = get.ClassNames(model)
            cl = labels(model.PrivClassNames);
        end
    end

    methods (Access=protected)
        %%% OBJECT CONSTRUCTOR %%%
        function model = MultinomialRegression()
            % MultinomialRegression constructor. Not intended to be called
            % directly. Use FITMNR to fit a MultinomialRegression object
        end

        %%% FITTING METHODS INHERITED FROM PARAMETRICREGRESSION %%%
        function model = assignData(model, X, Y, nvArgs)
            % ASSIGNDATA Assign data for fitting
            %   This method assigns training data to object properties, and
            %   populates relevant metadata for it. It takes in:
            %       model - The MultinomialRegression object to set
            %               training data properties for
            %       X - Matrix or table of predictor data. If a table, it
            %           may also contain the weights and response
            %       Y - Response vector, matrix of class counts, string of
            %           a formula, or the name of a column in X to treat as
            %           the response
            %       nvArgs - A struct of the name-value arguments for the
            %                fit process
            w = nvArgs.Weights;

            [yOrig, yNumer] = MultinomialRegression.convertYToNumeric(X, Y, model.Formula.ResponseName);
            if istable(X)
                if ~internal.stats.isString(Y)
                    % Given table X and response vector Y
                    % Add Y to the X table
                    X = [X, table(yNumer, 'VariableNames', {model.Formula.ResponseName})];
                else
                    % Given table X and formula/column name Y
                    % Replace the table entry with the numeric version
                    X.(model.Formula.ResponseName) = yNumer;
                end
            end

            asCat = nvArgs.CategoricalPredictors;
            if islogical(asCat) && isvector(asCat)
                % Logical vector
                if isfloat(X) && (numel(asCat) ~= model.Formula.NPredictors)
                    % Logical vector must be the same as the number of
                    % columns of X
                    error(message('stats:MultinomialRegression:BadCategoricalPredictorSize', ...
                        model.Formula.NPredictors));
                elseif istable(X)
                    % The object expects the vector to be the same length
                    % as the number of predictors, but FitObject expects
                    % it to be the same as the number of variables total.
                    asCat = asCat(:);
                    if numel(asCat) == model.Formula.NPredictors
                        % Add extra logicals to pad out the vector
                        origAsCat = asCat;
                        asCat = false(1, model.Formula.NVars);
                        asCat(model.Formula.InModel) = origAsCat;

                        % The table may have variables in it that aren't
                        % part of the model - make sure to encode them as
                        % categorial or not properly
                        asCat(~model.Formula.InModel) = varfun(@internal.stats.isDiscreteVar,...
                            X(:,~model.Formula.InModel),'OutputFormat','uniform');
                    else
                        % Logical vector length isn't the same as the
                        % number of predictors or the number of variables
                        error(message('stats:MultinomialRegression:BadCategoricalPredictorSize', ...
                            model.Formula.NPredictors))
                    end
                end
            elseif strcmp(asCat, 'all')
                if istable(X)
                    % Pad categorical var
                    asCat = true(1, model.Formula.NPredictors+1);
                else
                    asCat = true(1, model.Formula.NPredictors);
                end
            elseif nvArgs.wasSet.CategoricalPredictors && ...
                    internal.stats.isStrings(asCat) && any(~ismember(asCat, model.Formula.PredictorNames))
                % Provdided a name for CategoricalPredictors that doesn't
                % appear as a predictor
                error(message('stats:classreg:learning:internal:utils:CategoricalNotPredictor'))
            end

            model = assignData@classreg.regr.ParametricRegression(model, X, yNumer, w, asCat, model.Formula.VariableNames, []);

            if isfloat(X)
                model.PredictorTypes = 'numeric';
            else
                model.PredictorTypes = 'mixed';
            end
            model.ResponseType = 'categorical';

            % assignData higher up can't handle some forms of responses that
            % this object can. Hold on to the original form of Y so that we
            % can update the object later
            model.WorkingValues.yOrig = yOrig;

            % Set some of the nvArgs that are relevant for fitting but
            % don't have properties to hold them on the object
            model.WorkingValues.EstDisp = nvArgs.EstimateDispersion;
            model.WorkingValues.IterLimit = nvArgs.IterationLimit;
            model.WorkingValues.Tol = nvArgs.Tolerance;
        end

        function model = selectObservations(model, exclude)
            % SELECTOBSERVATIONS Determine which observations to use in the fit
            %   SELECTOBSERVATIONS determines which observations should be
            %   used in the fit. Any observations that have missing data
            %   are excluded from the fit. It takes in:
            %       model - The model to fit
            %       exclude - A logical vector of observations to exclude.
            model = selectObservations@classreg.regr.ParametricRegression(model, exclude);
            subset = model.ObservationInfo.Subset;

            % Object still operating on numeric conversion of response
            % data. Update the original response data based on missing data
            yOrig = model.WorkingValues.yOrig;
            if size(yOrig,2) > 1
                % Matrix of class counts
                % If the matrix has NaN values, we need to include them as
                % part of the set of missing values
                missing = model.ObservationInfo.Missing | any(isnan(yOrig),2);
                subset = subset & ~missing;
                model.ObservationInfo.Missing = missing;
                model.ObservationInfo.Subset = subset;
            end
        end

        function model = fitter(model)
            % FITTER Fit a model to data
            %   FITTER actually performs the fit of the input model to a
            %   dataset. It takes in:
            %       model - The model to fit.
            %   This method is called in succession after other methods do
            %   assignment of training data to the object.
            X = getData(model);
            Y = getResponse(model);

            if size(model.WorkingValues.yOrig, 2) > 1
                % Given Y as a matrix of counts - use that instead
                Y = model.WorkingValues.yOrig;
            end
            [designMat, predNames] = designMatrix(model, X);

            % Fit the model, set some object properties
            [coefs, model.Deviance, stats] = mnrfit(designMat, Y, 'Link', model.Formula.Link, ...
                'Model', model.ModelType, 'Interactions', model.IncludeClassInteractions, ...
                'EstDisp', model.WorkingValues.EstDisp, 'IterationLimit', model.WorkingValues.IterLimit,...
                'Tolerance', model.WorkingValues.Tol, 'Weights', model.ObservationInfo.Weights); %#ok<MNRFIT>
            model.CoefSize = size(coefs);
            model.Coefs = coefs(:);
            model.CoefficientCovariance = stats.covb;
            model.DFE = stats.dfe;
            model.PrivFitStats = stats;
            model.WorkingValues.PredNames = predNames;
        end

        function model = postFit(model)
            % POSTFIT Perform post-fitting actions on the object
            %   POSTFIT updates the object to set properties after the
            %   object has been fitted to data. It takes in:
            %       model - A fitted model to perform post-fit actions on


            % Now that the fitting is complete, convert the numeric
            % representation of Y on the object back to its original type
            yNumer = getResponse(model);
            yOrig = model.WorkingValues.yOrig;

            % Update Variables property
            if fitFromDataset(model)
                model.Data.(model.ResponseName) = yOrig;
            else
                model.Data.y = yOrig;
            end

            % Update VariableInfo
            model.VariableInfo{model.RespLoc, 'IsCategorical'} = true;
            model.VariableInfo{model.RespLoc, 'Class'} = {class(yOrig)};

            if size(yOrig, 2) > 1
                % Matrix of counts - 1 class per column
                range = 1:size(yOrig, 2);
                yNumer = yOrig;
            else
                range = unique(yOrig)';
                range = range(~ismissing(range));
            end
            model.VariableInfo{model.RespLoc, 'Range'} = {range};

            % Set Some other properties
            model.PrivClassNames = classreg.learning.internal.ClassLabel(range);
            model.CoefficientNames = model.createCoefficientNames(model.WorkingValues.PredNames, ...
                model.IncludeClassInteractions, string(labels(model.PrivClassNames(:))), model.ModelType);
            model.LogLikelihood = model.getlogLikelihood();
            model.LogLikelihoodNull = model.logLikelihoodNull(yNumer);

            % Remove the Excluded vector - there is no way to specify
            % excluded data like in LM/GLM
            model.ObservationInfo.Excluded = [];
        end

        %%% INHERITED METHODS FROM PARAMETRICREGRESSION TO BE USED FOR OBJECT METHODS %%%
        function L0 = logLikelihoodNull(model, Y)
            if nargin == 1
                % Pull Y from the object
                Y = getResponse(model);
                if size(model.WorkingValues.yOrig, 2) > 1
                    % Given Y as a matrix of counts - use that instead
                    Y = model.WorkingValues.yOrig;
                else
                    % Convert Y to index
                    Y = grp2idx(Y);
                end
            end
            ws = warning('off', 'stats:mnrfit:NominalNoInteractions');
            restoreWarn = onCleanup(@() warning(ws));
            [~, dev] = mnrfit(zeros(size(Y,1),0), Y, 'Link', model.Link, ...
                'Model', model.ModelType, 'Weights', model.ObservationInfo.Weights, ...
                'Interactions', model.IncludeClassInteractions); %#ok<MNRFIT>
            L0 = -.5 * dev;
        end

        function out = get_diagnostics(~)
            % Used to populate a dynamic 'Diagnostics' property. We don't
            % support this on the object. Silently just give no answer,
            % instead of erroring
            out = [];
        end

        function out = getlogLikelihood(model)
            % GETLOGLIKELIHOOD Compute the log-likelihood of the model
            out = -.5 * model.Deviance;
        end

        function resid = get_residuals(model, type)
            % GETRESIDUALS Get the residuals of the fitted model
            %   GETRESIDUALS returns the residuals of a fitted
            %   MultinomialRegression object model. It takes in:
            %       model - A fitted MultinomialRegression object
            %       type - One of 'raw', 'pearson', or 'deviance',
            %              indicating which residual to return
            %   If type is not given, a table with all residual types is
            %   returned
            if nargin == 1
                % Generate residual table
                resid = table(model.PrivFitStats.resid, model.PrivFitStats.residp, ...
                    model.PrivFitStats.residd, 'VariableNames', {'Raw', 'Pearson', 'Deviance'});
                return
            end
            if strcmpi(type, 'raw')
                resid = model.PrivFitStats.resid;
            elseif strcmpi(type, 'pearson')
                resid = model.PrivFitStats.residp;
            elseif strcmpi(type, 'deviance')
                resid = model.PrivFitStats.residd;
            else
                error(message('stats:GeneralizedLinearModel:UnrecognizedResidual', type));
            end
        end

        function r2 = get_rsquared(model, type)
            %GET_RSQUARED Compute the ordinary and adjusted R^2 values
            %   GET_RSQUARED computes pseudo-R^2 values suggested by
            %   McFadden. By default, both the adjusted and ordinary R^2
            %   are computed, but TYPE can be provided to select either
            %   one.
            loglFull = model.LogLikelihood;
            loglNull = model.LogLikelihoodNull;
            ordR2 = 1 - (loglFull/loglNull);
            adjR2 = 1 - ((loglFull - model.NumCoefficients)/loglNull);
            if nargin < 2
                r2 = struct('Ordinary', ordR2, 'Adjusted', adjR2);
            elseif strcmp(type, 'Ordinary')
                r2 = ordR2;
            elseif strcmp(type, 'Adjusted')
                r2 = adjR2;
            else
                error(message('stats:MultinomialRegression:UnknownRSquared'));
            end
        end

        function tbl = tstats(model)
            % TSTATS Generate table of info about coefficients
            %   TSTATS creates a table with the fitted coefficients, their
            %   SE, t stats, and p values. It takes in:
            %       model - A fitted MultinomialRegression object to
            %               generate a table for
            coefs = model.Coefs(:);
            se = sqrt(diag(model.CoefficientCovariance));
            tbl = table(coefs, se, model.PrivFitStats.t(:), model.PrivFitStats.p(:), ...
                'VariableNames',{'Value' 'SE' 'tStat', 'pValue'}, ...
                'RowNames',model.CoefficientNames);
        end

        function varargout = predictPredictorMatrix(model, XPred)
            % PREDICTPREDICTORMATRIX Predict when given the design matrix
            %   XPred is a design matrix. This method is used to bypass
            %   prediction method processing when a table isn't given
            [varargout{1:nargout}] = model.mnrvalWrapper(XPred);
        end

        %%% OTHER HELPER METHODS FOR OBJECT METHODS %%%
        function [design, predNames] = designMatrix(model, X)
            % DESIGNMATRIX Construct a design matrix
            %   DESIGNMATRIX creates a design matrix for a set of
            %   predictors X, and returns it along with the predictor names.
            %       It takes in:
            %       model - MultinomialRegression object that informs the
            %               matrix to create
            %       X - Matrix/table of predictors to convert to a design
            %           matrix

            % Get the terms from the formula, but remove the intercept row.
            % mnrfit automatically includes an intercept row.
            terms = model.Formula.Terms;
            terms(all(terms==0,2),:) = [];
            if istable(X)
                predNames = X.Properties.VariableNames;
                if ~all(ismember(model.PredictorNames, predNames))
                    error(message('stats:classreg:regr:TermsRegression:MissingVariable'));
                end
                [tf,varLocs] = ismember(X.Properties.VariableNames,model.VariableNames);
                if ~all(tf)
                    % Remove extra variables from X so that we don't have to
                    % construct Terms, IsCategorical, Range, and
                    % DummyVarCoding values.
                    X = X(:,tf);
                end
                varLocs = varLocs(tf);
            else
                nvars = model.NumVariables;
                if size(X,2) ~= nvars-1
                    error(message('stats:classreg:regr:TermsRegression:WrongXColumns', nvars - 1));
                end
                varLocs = true(1,nvars);
                varLocs(model.RespLoc) = false;
            end
            [design,~,~,~,predNames] = classreg.regr.modelutils.designmatrix(X,'Model',terms(:,varLocs), ...
                'CategoricalVars',model.VariableInfo.IsCategorical(varLocs), ...
                'VarNames',model.Formula.VariableNames(varLocs), 'Intercept', false, ...
                'CategoricalLevels', model.VariableInfo.Range(varLocs));
        end

        function [out, probs, lower, upper] = mnrvalWrapper(model, designMat, varargin)
            % MNRVALWRAPPER A wrapper around the mnrval method
            %   MNRVALWRAPPER is used to call mnrval, which does prediction
            %   for a fitted model. It takes in:
            %       model - Fitted MultinomialRegression object
            %       designMat - The design matrix to predict on
            %       varargin - NV args supported by mnrval
            coefs = reshape(model.Coefs, model.CoefSize);
            args = {'Type', 'Confidence'};
            defs = {'category', .05};
            [type, conf] = internal.stats.parseArgs(args, defs, varargin{:});

            computeCI = (nargout > 2);
            if strcmpi(type, 'category')
                % Don't need to recompute probabilities - get them all at
                % once with CI if requested
                if computeCI
                    [catProbs, lower, upper] = mnrval(coefs, designMat, model.PrivFitStats, ...
                        'Model', model.ModelType, 'interactions', model.IncludeClassInteractions, ...
                        'link', model.Formula.Link, 'confidence', conf); %#ok<MNRVAL>
                else
                    catProbs = mnrval(coefs, designMat, 'Model', model.ModelType, ...
                        'interactions', model.IncludeClassInteractions, 'link', model.Formula.Link); %#ok<MNRVAL>
                end
                probs = catProbs;
            else
                % Need to get categorical probabilities first to predict
                % the class label, then recompute the
                % conditional/cumulative probabilities
                catProbs = mnrval(coefs, designMat, 'Model', model.ModelType, ...
                    'interactions', model.IncludeClassInteractions, 'link', model.Formula.Link); %#ok<MNRVAL>

                if computeCI
                    [probs, lower, upper] = mnrval(coefs, designMat, model.PrivFitStats, ...
                        'Model', model.ModelType, 'interactions', model.IncludeClassInteractions, ...
                        'link', model.Formula.Link, 'type', type, 'confidence', conf); %#ok<MNRVAL>
                else
                    probs = mnrval(coefs, designMat, model.PrivFitStats, ...
                        'Model', model.ModelType, 'interactions', model.IncludeClassInteractions, ...
                        'link', model.Formula.Link, 'type', type, 'confidence', conf); %#ok<MNRVAL>
                end
            end

            if computeCI
                % Put lower and upper into expected form
                lower = probs - lower;
                lower(lower < 0) = 0;
                upper = probs + upper;
                upper(upper > 1) = 1;
            end

            [~, inds] = max(catProbs, [], 2);
            out = model.ClassNames(inds);
        end

        function [scores, xVals] = predictSliceScores(model, variedPredName, preddata)
            % PREDICTSLICESCORES Predict scores for fixed predictor values
            %   PREDICTSLICESCORES Predicts the scores for use in
            %   generating slice plots. It fixes all predictors but one,
            %   and varies the remaining one to generate a set of scores.
            %   It takes in:
            %       model - The model to predict with
            %       variedPredName - The name of the predictor to vary to
            %                        generate the score surface
            %       preddata - Predictor data used to set the limits of the axes. 
            %                  If preddata is empty, the limits are the range 
            %                  of the training data.
            % Determine the number of datapoints we need in the plot
            varInfo = model.VariableInfo(variedPredName, :);
            variedPredName = cellstr(variedPredName);
            ranges = classreg.regr.modelutils.getDataRangesForPredictorData(model, preddata);
            [~, rangeInd] = ismember(variedPredName, model.PredictorNames);
            if istable(preddata)
                % Remove any extra columns and align the table order
                preddata = preddata(:, model.PredictorNames);
            end

            varRange = ranges{rangeInd};
            if varInfo.IsCategorical
                % For a categorical predictor, we show all categories in
                % preddata
                xVals = varRange;
                if ischar(xVals)
                    xVals = cellstr(xVals);
                end
                numPoints = numel(xVals);
            else
                % For a continuous predictor, we show a uniform grid of 40
                % points
                numPoints = 40;
                xVals = linspace(varRange(1), varRange(2), numPoints)';
            end
            data = table(xVals(:), 'VariableNames', variedPredName);

            % Iterate over all the predictors to generate a table with
            % vectors of average/baseline predictor values
            havePreddata = ~isempty(preddata);
            for i = 1:model.NumPredictors
                currPred = model.PredictorNames(i);
                if ~strcmp(currPred, variedPredName)
                    % Table already has varied predictor in it
                    if model.VariableInfo.IsCategorical(currPred)
                        % Dummy variable coding is 'reference' - use the first
                        % categorical level
                        if ischar(ranges{i})
                            val = ranges{i}(1,:);
                        else
                            val = ranges{i}(1);
                        end
                    else
                        % Continuous predictor - use the mean
                        if havePreddata
                            % Get the mean of the training data
                            if istable(preddata)
                                val = mean(preddata{:, currPred});
                            else
                                val = mean(preddata(:,i));
                            end
                        else
                            val = mean(model.Variables{model.ObservationInfo.Subset, currPred});
                        end
                    end
                    data = [data, table(repmat(val, numPoints, 1), 'VariableNames', currPred)]; %#ok<AGROW>
                end
            end
            [~, scores] = model.predict(data);

            if any(isnan(scores), 'all')
                % Given bad data, most likely a category that doesn't exist
                error(message('stats:classreg:regr:modelutils:BadPredDataPredictNaN'))
            end
        end
    end

    methods (Hidden, Static)
        %%% OBJECT FITTING HELPER FUNCTIONS %%%
        function model = fit(X, Y, varargin)
            % FIT fit a MultinomialRegression object. Not intended to be called
            % directly. Use FITMNR to fit a MultinomialRegression
            %
            %   See also fitmnr, fitglm

            % Validate and extract inputs
            MultinomialRegression.validateX(X);
            Y = MultinomialRegression.validateY(Y, istable(X), size(X, 1));
            nvArgs = MultinomialRegression.parseAndValidateNVArgs(varargin);

            % Create model
            model = MultinomialRegression();
            model.IncludeClassInteractions = nvArgs.IncludeClassInteractions;
            model.ModelType = nvArgs.ModelType;
            model.Formula = MultinomialRegression.createFormula(X, Y, nvArgs);

            % Assign data and fit the model, using methods in the
            % ParametricRegression hierarchy
            model = assignData(model, X, Y, nvArgs);
            model = doFit(model);
        end

        function fitFcnName = getFitFcnName()
            % Return the name of the public function to fit a
            % MultinomialRegression object
            fitFcnName = 'fitmnr';
        end

        function [yOrig, yNumer] = convertYToNumeric(X, Y, respName)
            % CONVERTYTONUMERIC Convert a response to a numeric vector
            %   CONVERTYTONUMERIC takes a response vector, the name of
            %   a column in X, or a formula Y, and converts it to numeric.
            %   It also returns the original response. It takes in:
            %       X - The predictor data to fit. Either a double/single
            %           matrix or table
            %       Y - The response vector. If X is a table, Y can also be
            %           a formula in Wilkinson notation, or the name of a
            %           column in X
            %       respName - The name of the response vector in fitting
            %       the object

            % Convert from categorical Y to numeric Y
            if internal.stats.isString(Y)
                % X must be a table, Y is in the table
                % Pull it out and validate its type.
                yOrig = X.(respName);
                yOrig = MultinomialRegression.validateY(yOrig, false, size(X, 1));
            else
                % Y is a separate response vector
                yOrig = Y;
            end
            if isfloat(yOrig) && (size(yOrig, 2) > 1)
                % Matrix of counts
                % Use the first column as a placeholder for yNumer, as that
                % is sufficient enough for some basic validation
                yNumer = yOrig(:,1);
            else
                yNumer = grp2idx(yOrig);
            end
        end

        function formula = createFormula(X, Y, nvArgs)
            % CREATEFORMULA Create a formula that represents the model
            %   CREATEFORMULA creates a classreg.regr.LinearFormula
            %   representing the fitted MultinomialRegression object. It
            %   takes in:
            %       X - The predictor data. A double/single matrix or table
            %       Y - The response data. If X is a table, can also be a
            %           string of a column in X, or a formula in Wilkinson
            %           notation
            %       nvArgs - A struct of NV args used to fit the
            %                MultinomialRegression object

            % Before creating the formula, we need to determine what the
            % names of the variables in the universe are
            modelSpec = 'linear';
            if istable(X)
                % Y can be:
                % 1. A response vector/matrix
                % 2. A name of a column in X
                % 3. A formula
                varNames = X.Properties.VariableNames;
                if internal.stats.isString(Y)
                    % Y is a formula or the name of a table column
                    if ismember(Y, varNames)
                        % Column name
                        responseName = Y;
                        if nvArgs.wasSet.PredictorNames
                            modelSpec = {'linear', nvArgs.PredictorNames};
                        end
                    else
                        % Formula - presumably
                        if ~contains(Y, '~')
                            % Neither a valid formula nor a table column
                            % name. Assume that we were given a bad table
                            % column name
                            error(message('stats:MultinomialRegression:NoMatchingResponseName', Y));
                        end
                        responseName = [];
                        modelSpec = Y;
                    end
                else
                    % Response Y - pull response name from nvArgs, or create
                    % a default name
                    if nvArgs.wasSet.ResponseName
                        responseName = nvArgs.ResponseName;
                    else
                        responseName = matlab.lang.makeUniqueStrings([varNames, 'y']);
                        responseName = responseName{end};
                    end
                    varNames = [varNames, responseName];
                    if nvArgs.wasSet.PredictorNames
                        % Given a subset of table variables to use
                        modelSpec = {'linear', nvArgs.PredictorNames};
                    end
                end
            else
                % Matrix X - either we're given a set of names, or we
                % generate a default set
                [varNames, ~, responseName] = ...
                    classreg.regr.FitObject.getVarNames([], nvArgs.PredictorNames, nvArgs.ResponseName, size(X, 2));
            end
            formula = classreg.regr.LinearFormula(modelSpec, varNames, responseName, true);

            if strcmpi(formula.Link, 'identity')
                % No link specified in formula - use the NV arg for it
                formula.Link = nvArgs.Link;
            elseif nvArgs.wasSet.Link && ~strcmpi(formula.Link, nvArgs.Link)
                % User specified one kind of link in the formula, and a
                % different with the 'Link' NV pair
                error(message('stats:MultinomialRegression:InconsistentLinkFcn'));
            end
            validatestring(formula.Link, {'logit', 'probit', 'loglog', 'comploglog'}, '', 'Link');

            if nvArgs.wasSet.ResponseName
                if ~strcmp(nvArgs.ResponseName, formula.ResponseName)
                    % Ensure ResponseName NV arg aligns with the formula
                    % User can, e.g., provide one response name in a
                    % formula and another in the ResponseName NV arg
                    error(message('stats:classreg:learning:internal:utils:ConflictingArg', 'ResponseName'));
                end
            end
        end

        function coefNames = createCoefficientNames(predNames, inclInter, classNames, modelType)
            predNames = string(predNames(:));
            if inclInter
                % Coefficients are arranged with predictors on the rows and
                % classes on the columns. Predictors are different across
                % each class, and required extra naming to distinguish them
                coefNames = ["(Intercept_" + classNames(1:end-1)' + ")"; predNames + "_" + classNames(1:end-1)'];
            else
                % No interactions - coeffs are a single vector. The first
                % k-1 rows are for intercept terms, and the rest of the
                % entries are each of the predictors. Predictors are the
                % same across categories, and don't need additional naming
                if strcmpi(modelType, 'nominal')
                    % Nominal model without class interactions is just
                    % intercept terms
                    coefNames = "(Intercept_" + classNames(1:end-1) + ")";
                else
                    coefNames = ["(Intercept_" + classNames(1:end-1) + ")"; predNames];
                end
            end
            coefNames = coefNames(:);
        end

        %%% INPUT VALIDATION FUNCTIONS %%%
        function validateX(X)
            % VALIDATEX Validate predictor data
            %   VALIDATEX validates predictor data X. It takes in:
            %       X - A double/single matrix or table or predictor data

            % Do basic type validation. Further validation happens later
            if istall(X) || isa(X, 'gpuArray')
                error(message('stats:MultinomialRegression:TallGPUPredictorsNotSupported'));
            end
            if isfloat(X)
                validateattributes(X, {'double', 'single'}, {'real', '2d', 'nonsparse'}, '', 'X');
            elseif ~istable(X)
                error(message('stats:MultinomialRegression:BadPredictorType'));
            end
        end

        function Y = validateY(Y, xIsTable, numXRows)
            % VALIDATEY Validate response data
            %   VALIDATEY Performs basic type validation on Y, a
            %   representation of the response variable. It takes in:
            %       Y - Can be:
            %               1. Label type - categorical, cellstr/string, numeric, logical
            %               2. Matrix of class counts - double/single/logical
            %               3. Char vector/string of a column name/formula -
            %                  only supported if X is a table
            %       xIsTable - Logical indicator if predictor data is in a table
            %       numXRows - Number of rows of predictors
            validateParams = {'', 'Y'};
            % Treat logical as numeric
            if isfloat(Y) || islogical(Y)
                % Either (1) or (2)
                if isvector(Y)
                    % (1)
                    validateattributes(Y, {'double', 'single', 'logical'}, {'real', 'nonsparse'}, validateParams{:})
                else
                    % (2)
                    % Only validate rows without NaN data
                    validateattributes(Y, {'double', 'single', 'logical'}, {'2d'}, validateParams{:});
                    nonNan = Y(all(~isnan(Y), 2), :);
                    validateattributes(nonNan, {'double', 'single', 'logical'}, {'real', 'nonsparse', 'integer', 'nonnegative'}, validateParams{:});
                end
            elseif internal.stats.isString(Y)
                % (3), or (1) but there is only one row of matrix predictor data
                if (numXRows == 1)
                    Y = cellstr(Y);
                elseif ~xIsTable
                    error(message('stats:MultinomialRegression:BadYType'))
                end
            else
                if ischar(Y) || isstring(Y)
                    Y = cellstr(Y);
                end
                if ~((iscellstr(Y) || iscategorical(Y)) && isvector(Y))
                    error(message('stats:MultinomialRegression:BadYType'))
                end
            end
        end

        function nvArgs = parseAndValidateNVArgs(argList)
            % PARSEANDVALIDATENVARGS Parse and validate fit NV args
            %   PARSEANDVALIDATENVARGS parses the set of NV args supported
            %   in fitting the object, and performs basic validation on
            %   some. It takes in:
            %       argList - A cellstr of NV args
            [argList{:}] = convertStringsToChars(argList{:});
            args = {'Link', 'ModelType', 'IncludeClassInteractions', 'EstimateDispersion', ...
                'Tolerance', 'IterationLimit', 'Weights', 'CategoricalPredictors', ...
                'ResponseName', 'PredictorNames'};
            defs = {'logit', 'nominal', [], false, 1e-6, 100, [], [], [], []};
            [link, modelType, incInter, estDisp, tol, numIter, wts, catPreds, respName, predNames, wasSet] = ...
                internal.stats.parseArgs(args, defs, argList{:});

            % Fully validate the simpler NV args
            fitFcn = MultinomialRegression.getFitFcnName();
            link = validatestring(link, {'logit', 'probit', 'comploglog', 'loglog'}, fitFcn, 'Link');
            modelType = validatestring(modelType, {'nominal', 'ordinal', 'hierarchical'}, fitFcn, 'ModelType');
            if ~wasSet.IncludeClassInteractions
                if strcmpi(modelType, 'ordinal')
                    incInter = false;
                else
                    incInter = true;
                end
            else
                incInter = internal.stats.parseOnOff(incInter, 'IncludeClassInteractions');
            end
            estDisp = internal.stats.parseOnOff(estDisp, 'EstimateDispersion');
            validateattributes(tol, {'double', 'single'}, {'scalar', 'real', 'positive', 'nonsparse'}, ...
                fitFcn, 'Tolerance');
            validateattributes(numIter, {'double', 'single'}, {'real', 'positive', 'integer', 'nonsparse'}, ...
                fitFcn, 'IterationLimit');

            % Partially validate PredictorNames and ResponseName. More
            % complex validation is done later in fitting
            if wasSet.PredictorNames
                if ischar(predNames)
                    predNames = cellstr(predNames);
                end
                if ~iscellstr(predNames) || ~isvector(predNames)
                    error(message('stats:classreg:learning:FullClassificationRegressionModel:prepareDataCR:BadPredictorType'));
                end
            end

            if wasSet.ResponseName
                validateattributes(respName, {'char', 'string'}, {'scalartext'}, fitFcn, 'ResponseName');
            end

            if wasSet.Weights
                validateattributes(wts, {'double', 'single'}, {'vector', 'nonnegative', 'real', 'nonsparse'},'fitmnr', 'Weights');
            end

            % Weights and CategoricalPredictors are fully validated later
            % in fitting
            nvArgs = struct();
            nvArgs.Link = link;
            nvArgs.ModelType = modelType;
            nvArgs.IncludeClassInteractions = incInter;
            nvArgs.EstimateDispersion = estDisp;
            nvArgs.Tolerance = tol;
            nvArgs.IterationLimit = numIter;
            nvArgs.Weights = wts;
            nvArgs.CategoricalPredictors = catPreds;
            nvArgs.ResponseName = respName;
            nvArgs.PredictorNames = predNames;
            nvArgs.wasSet = wasSet;
        end

        function validateClassName(className, allowedSize, paramName)
            % VALIDATECLASSNAME Validates a class name
            %   VALIDATECLASSNAME validates that a given class name is the
            %   proper type and size. It takes in:
            %       - className - The class name. It must be a label type
            %       - allowedSize - 'scalar' or 'vector', indicating if
            %                        multiple class names can be given, or
            %                        just one
            %      - paramName -  The name of the parameter that the class
            %                     name corresponds to
            if ischar(className)
                className = cellstr(className);
            end
            if strcmpi(allowedSize, 'vector')
                sizeFcn = @(x) isvector(x);
            else
                sizeFcn = @(x) isscalar(x);
            end
            if ~(iscellstr(className) && sizeFcn(className))
                validateattributes(className, {'logical', 'categorical', 'string', 'double', 'single'}, ...
                    {'nonempty', allowedSize}, '', paramName);
            end
        end

        function [numTrials, type] = parseAndValidateRandomArgs(varargin)
            % parseAndValidateRandomArgs Parse and validate extra args to
            % random
            %   parseAndValidateRandomArgs parses and validates the extra
            %   arguments given to the MultinomialRegression.random method.
            %   There are two possible arguments, which can appear in any
            %   order (or may not appear at all):
            %       - numTrials - The number of trials to generate random
            %                     samples for. Default is 1.
            %       - type - The type that the random sample output should
            %                appear in. The options are 'counts' or
            %                'samples' (default).
            if nargin == 0
                % None given - use defaults
                numTrials = 1;
                type = "samples";
            elseif nargin == 1
                % Either numTrials or type is given
                if isfloat(varargin{1})
                    % Assume numTrials
                    numTrials = varargin{1};
                    type = "samples";
                else
                    numTrials = 1;
                    type = varargin{1};
                end
            elseif nargin == 2
                % Both numTrials and type given
                if isfloat(varargin{1})
                    numTrials = varargin{1};
                    type = varargin{2};
                else
                    numTrials = varargin{2};
                    type = varargin{1};
                end
            else
                % Too many args
                error(message('MATLAB:maxrhs'));
            end
            validateattributes(numTrials, {'double', 'single'}, ...
                {'real', 'integer', 'scalar', 'nonsparse', 'positive'}, '', 'numTrials');
            type = validatestring(type, ["counts", "samples"], '', 'type');
        end

        %%% PLOT METHOD HELPERS %%%
        function [parent, model, args] = parseParentAndModel(varargin)
            % PARSEPARENTANDMODEL Parses out the axes, object, and NV args
            %   PARSEPARENTANDMODEL parses an optional axes/figure, required 
            %   object and any NV args from the input varargin. The input 
            %   is a cell array that at minimum contains the object. It 
            %   may also contain an axes/figure object as the first input
            [varargin{:}] = convertStringsToChars(varargin{:});
            [parent, varargin] = axescheck(varargin{:});
            if isa(varargin{1}, 'matlab.ui.Figure')
                parent = varargin{1};
                model = varargin{2};
                args = varargin(3:end);
            else
                model = varargin{1};
                args = varargin(2:end);
            end
        end

        function f = plotInteractive(model, preddata, parent, NVArgs)
            % PLOTINTERACTIVE Generate an interactive slice plot
            %   PLOTINTERACTIVE creates an interactive plot of the slices
            %   of the regression surface defined by model. It takes in:
            %       model - The MultinomialRegression object to generate slices for
            %       preddata - Predictor data to use to generate the
            %                  slices. May be empty.
            %       parent - A parent figure to use when creating the
            %                plotSlice graphic
            %       NVArgs.ClassToPlot - An NV arg determine which class to
            %                            generate the slices for. Valid
            %                            choices are model.ClassNames
            arguments
                model(1,1) MultinomialRegression
                preddata = []
                parent = []
                NVArgs.ClassToPlot {MultinomialRegression.validateClassName(NVArgs.ClassToPlot, 'scalar', 'ClassToPlot')} = model.ClassNames(1);
            end
            cl = classreg.learning.internal.ClassLabel(NVArgs.ClassToPlot);
            [~, ind] = ismember(cl, model.PrivClassNames);
            if ind == 0
                error(message('stats:MultinomialRegression:UnrecognizedClass', char(cl)));
            end
            if ~isempty(parent) && ~isa(parent, 'matlab.ui.Figure')
                % Given an axes handle
                error(message('stats:classreg:regr:modelutils:IncorrectFigureType'));
            end
            f = classreg.regr.modelutils.plotSlice(model, 'ClassInd', ind, ...
                'IsClassif', true, 'PredictorData', preddata, 'HasSimulCI', false, ...
                'Parent', parent);
        end

        function hOut = plotHistogram(model, ax, plottype, preddata, varargin)
            % PLOTHISTOGRAM Generate a slice plot represented by a histogram
            %   PLOTHISTOGRAM creates a slice plot of the regression
            %   surface in the fitted MultinomialRegression object model.
            %   It can display this surface in two styles. It takes in:
            %       model - The MultinomialRegression object to plot
            %       ax - An Axes on which to plot. May be empty
            %       plottype - 'histogram' or 'stackedhist', which
            %                   generates a histogram of slices per class
            %                   overlaid on eachother or stacked on eachother,
            %                   respectively
            %       preddata - A set of data to generate predictions from.
            %                  This is used to set the axes limits and the
            %                  reference values when generating the plot.
            %                  It may be empty, in which case the ranges
            %                  from the training data are used
            %       varargin - A set of NV args. Arguments to set graphics
            %                  properties for Histogram or Bar objects are
            %                  supported.
            %                  A named argument, 'PredictorToVary' is also
            %                  supported. It determines which predictor to
            %                  vary on the x axis. It must be one of
            %                  model.PredictorNames
            [predName, ~, plotArgs] = internal.stats.parseArgs({'PredictorToVary'}, {model.PredictorNames(1)}, varargin{:});
            if ~internal.stats.isString(predName,true)
                error(message('stats:MultinomialRegression:BadPredictorNameSpecified'));
            elseif ~ismember(predName, model.PredictorNames)
                error(message('stats:MultinomialRegression:NoMatchingPredictorName', predName))
            end
            [scores, xVals] = model.predictSliceScores(predName, preddata);

            if isempty(ax)
                ax = gca();
            end
            if strcmp(plottype, 'stackedhist')
                % Stacked histogram
                hOut = bar(ax, xVals, scores, 'stacked', plotArgs{:});
                ylim(ax, [0, 1.01]);
            else
                % Overlaid histogram
                nextPlot = ax.NextPlot;
                numClasses = size(scores, 2);
                hOut = gobjects(1, numClasses);
                isCat = ~isnumeric(xVals);

                % Add an extra point to xVals for numeric predictors
                % The way that histogram is called requires us to supply
                % the bin edges and the values of the bins. There must be 1
                % more edge than values of the bins
                if ~isCat
                    xVals = [xVals', xVals(end) + (xVals(end) - xVals(end-1))];
                end

                for i = 1:size(scores, 2)
                    if isCat
                        hOut(i) = histogram(ax, 'Categories', xVals, 'BinCounts', scores(:,i)', plotArgs{:});
                    else
                        hOut(i) = histogram(ax, 'BinEdges', xVals, 'BinCounts', scores(:,i)', plotArgs{:});
                    end
                    hold(ax, 'on');
                end
                ax.NextPlot = nextPlot;
            end
            xlabel(ax, predName, 'Interpreter', 'none');
            ylabel(ax, getString(message('stats:MultinomialRegression:label_Probability')));
            classNames = cellstr(model.PrivClassNames);
            for i = 1:numel(hOut)
                hOut(i).DisplayName = classNames{i};
            end
        end
    end

    methods
        function [p, stat] = testDeviance(MnrModel)
            %TESTDEVIANCE Analysis of deviance
            %   P = TESTDEVIANCE(MNR) returns the p value for an analysis of
            %   deviance table for the MultinomialRegression MNR. The p value
            %   is for a test of whether the fitted model fits significantly
            %   better than a constant model.
            %
            %   [P, STAT] = TESTDEVIANCE(MNR) returns the test statistic
            %   STAT used to generate the P value P.
            %
            %   Example: Fit a regression model and examine the significance
            %            of the fit
            %   load fisheriris
            %   mnr = fitmnr(meas,species);
            %   [p, stat] = testDeviance(mnr)
            %
            %   See also coefCI, coefTest
            arguments
                MnrModel(1,1) MultinomialRegression
            end
            devNull = -2*MnrModel.LogLikelihoodNull;
            devFull = MnrModel.Deviance;
            % DF - difference in number of parameters in full model vs
            % intercept model
            dfDiff = MnrModel.NumCoefficients - numel(MnrModel.ClassNames) + 1;

            if MnrModel.EstimateDispersion
                stat = max(0, devNull - devFull) / (MnrModel.Dispersion * dfDiff); % F statistic
                p = fcdf(stat, dfDiff, MnrModel.DFE, 'upper'); % fpval(fstat,dfr,dfe);
            else
                stat = max(0, devNull - devFull); % chi-squared statistic
                p = chi2cdf(stat, dfDiff, 'upper');
            end
        end

        function ysim = random(model, X, varargin)
            %RANDOM Generate random response values given predictor values.
            %   YSIM = RANDOM(MNR,X) generates a vector YSIM of random values
            %   from the MultinomialRegression MNR using predictor variables
            %   from the matrix/table X. X must contain all of the predictor
            %   variables used to create MNR. The output YSIM is computed by
            %   creating predicted values and generating response values from
            %   the fitted distribution.
            %
            %   YSIM = RANDOM(__, NUMTRIALS) generates a matrix YSIM
            %   with NUMTRIALS columns. For each row of X, it simulates
            %   generating NUMTRIALS pulls of response values from a multinomial
            %   distribution.
            %
            %   YSIM = RANDOM(__, TYPE) generates random data YSIM based on
            %   TYPE. TYPE can be 'samples' (default) or 'counts'. If
            %   'samples', YSIM is a matrix of random response values. If
            %   TYPE is 'counts', YSIM is a double matrix with size(X, 1)
            %   rows and numel(MNR.ClassNames) columns. YSIM(i, j) represents
            %   the number of times MNR.ClassNames(j) was pulled for observation
            %   X(i,:).
            %
            %   Example:
            %   % Fit the model to the iris dataset, but modify the
            %   % predictor data to make the problem more complex
            %   load fisheriris
            %   meas = meas(:, 1:2) + rand(150, 2);
            %
            %   % Fit the model, and generate 5 random samples for 3
            %   observations, one each of setosa, versicolor, and virginica.
            %   mnr = fitmnr(meas, species);
            %   numTrials = 5;
            %   samps = random(mnr, meas([5, 55, 105], :), numTrials)
            %
            %   See also MultinomialRegression, predict.
            model.validateX(X)
            [numTrials, type] = model.parseAndValidateRandomArgs(varargin{:});
            [~, probs] = model.predict(X);
            pulls = mnrnd(numTrials, probs);
            if strcmpi(type, "counts")
                % pulls is in counts format - just return that instead
                ysim = array2table(pulls, 'VariableNames', string(labels(model.PrivClassNames)));
            else
                % convert from counts to specific samples
                if islogical(model.ClassNames)
                    ysim = true(0,0);
                else
                    ysim = [];
                end
                elem = 1:numel(model.PrivClassNames);
                for i = 1:size(pulls, 1)
                    ysim = [ysim; model.ClassNames(repelem(elem, pulls(i,:)))']; %#ok<AGROW>
                end
            end
        end

        function varargout = predict(model, X, NVArgs)
            %PREDICT Predict response of MultinomialRegression
            %   YPRED = PREDICT(MNR, X) returns predicted class labels
            %   YPRED for a MultinomialRegression MNR and predictors X. X
            %   can be a matrix or table of predictor data. It must contain
            %   all predictors used to train the MNR object. YPRED has the
            %   same type as MNR.ClassNames.
            %
            %   [YPRED, PROBS] = PREDICT(MNR, X) returns category probabilities
            %   PROBS. PROBS is an N x K matrix numeric matrix for N observations
            %   and K classes.
            %
            %   [YPRED, PROBS, LOWER, UPPER] = PREDICT(MNR, X) returns
            %   the lower and upper confidence bounds in LOWER and UPPER,
            %   respectively, for the PROBS. LOWER and UPPER are N x K
            %   numeric matrices of 95% confidence intervals for each
            %   element of PROBS.
            %
            %   PREDICT(MNR, X, Name1=Val1, ...) specifies optional
            %   name-value arguments:
            %       'Alpha' - A double/single scalar the specifies the
            %                 confidence level of the confidence intervals
            %                 as 100*(1-ALPHA)%. The default is .05
            %
            %       'ProbabilityType' - One of 'category' (default),
            %                           'cumulative', or 'conditional',
            %                           indicating which kind of probability
            %                           to return in PROBS.
            %
            %   Example: Predict labels on training set
            %   load fisheriris
            %   mnr = fitmnr(meas, species);
            %   [labels, probs] = predict(mnr, meas)
            %
            %   See also random, partialDependence, coefCI, feval
            arguments
                model(1,1) MultinomialRegression
                X {MultinomialRegression.validateX(X)} = model.getData
                NVArgs.Alpha(1,1) {anova.validation.validateAlpha(NVArgs.Alpha, 'predict')} = .05
                NVArgs.ProbabilityType(1,1) string = "category"
            end
            prob = validatestring(NVArgs.ProbabilityType, ["category", "cumulative", "conditional"]);
            designMat = model.designMatrix(X);
            [varargout{1:nargout}] = model.mnrvalWrapper(designMat, 'type', prob, 'confidence', NVArgs.Alpha);
        end

        function [pd,varargout] = partialDependence(model,features,labels,varargin)
            %PARTIALDEPENDENCE Partial Dependence for one and two variables
            %   PD = partialDependence(MODEL,VARS,LABELS) computes the partial
            %   dependence PD between the predictor variables listed in VARS and the
            %   scores for the classes specified in LABELS by using the
            %   MultinomialRegression model MODEL and by averaging over the
            %   data used in training the model. You can specify one or two
            %   f
            %   variables using VARS. VARS is either column indices of predictors
            %   or names of predictors. The names must match the entries in
            %   MODEL.PredictorNames. You can specify one or multiple classes
            %   using LABELS. The values in LABELS must match and have the same
            %   data type as the class names in MODEL.ClassNames. If there
            %   are multiple values in both VARS and LABELS, then PD is a
            %   three-dimensional array.
            %
            %   PD = partialDependence(...,DATA) specifies the data to be used
            %   for averaging. DATA is a matrix or table of data to be used
            %   in place of the data used in fitting the model.
            %
            %   PD = partialDependence(...,Name=Value) specifies additional
            %   options using one or more name-value arguments:
            %
            %      'NumObservationsToSample':   An integer specifying the number
            %                                   of rows to sample at random from
            %                                   the dataset (either the data
            %                                   specified by argument DATA or the
            %                                   data used to train MODEL). The
            %                                   default is to use all rows.
            %
            %      'QueryPoints':               The points at which to calculate
            %                                   the partial dependence. For the
            %                                   case of a single predictor, the
            %                                   value of 'QueryPoints' must be a
            %                                   column vector. For the case of two
            %                                   predictors, the value of 'QueryPoints'
            %                                   must be a 1x2   cell array containing
            %                                   a separate vector for each predictor.
            %                                   The default is to use 100 points
            %                                   equally spaced across the range of
            %                                   the predictor.
            %
            %      'UseParallel'                Can be true to instruct the method to
            %                                   perform the averaging calculations in
            %                                   parallel (using parfor), or false
            %                                   (default) to instruct the method not
            %                                   to use parallel computation.
            %
            %   [PD,X,Y] = partialDependence(...) also returns the arrays X and Y
            %   containing the query point values of the first and second predictor
            %   in VARS, respectively.
            %
            %   Examples:
            %       load census1994
            %       X = adultdata(:,{'age','workClass','education_num','marital_status','race',...
            %           'sex','capital_gain','capital_loss','hours_per_week','salary'});
            %       f = fitmnr(X,'salary');
            %
            %       pd = partialDependence(f,'age','<=50K');
            %
            %       % Obtain values x of the variable 'age'
            %       [pd,x] = partialDependence(f,'age',{'<=50K','>50K'});
            %
            %       % plot partial dependence of the first and second classes scores on
            %       % 'age'
            %       plot(x,pd)
            %
            %       % Obtain values x and y of variables 'age' and 'workClass'
            %       [pd,x,y] = partialDependence(f,{'age','workClass'},{'<=50K','>50K'});
            %
            %       % plot partial dependence for class '<=50K' on variables
            %       % 'age' and 'workClass'
            %       surf(x,y,squeeze(pd(1,:,:)))
            %
            %       % plot partial dependence for class '>50K' on variables
            %       % 'age' and 'workClass'
            %       surf(x,y,squeeze(pd(2,:,:)))
            %
            %       % With additional Data
            %       pd = partialDependence(f,'age','<=50K',X);
            %       pd = partialDependence(f,1,'<=50K',X);
            %
            %       % With optional name-value pairs
            %       pd = partialDependence(f,1,'<=50K',X,'NumObservationsToSample',100);
            %       pd = partialDependence(f,1,'<=50K',X,'UseParallel',true);
            %
            %       % Provide alternative query points
            %       xi = linspace(min(X{:,'age'}),max(X{:,'age'}))';
            %       pd = partialDependence(f,1,'<=50K','QueryPoints',xi);
            %
            %       % Provide alternative query points
            %       xi = linspace(min(X{:,'age'}),max(X{:,'age'}))';
            %       pd = partialDependence(f,1,'<=50K','QueryPoints',xi);
            %
            %       xi = cell(1,2);
            %       xi{1} = linspace(min(X{:,1}),max(X{:,1}))';
            %       xi{2} = linspace(min(X{:,9}),max(X{:,9}))';
            %       pd = partialDependence(f,[1,9],'<=50K','QueryPoints',xi);
            %
            %   See also plotPartialDependence, plotSlice, predict
            features = convertStringsToChars(features);
            [varargin{:}] = convertStringsToChars(varargin{:});
            p = inputParser;
            addRequired(p,'Model');
            addRequired(p,'Var');
            addRequired(p,'Label');
            addOptional(p,'Data',getData(model)); % Default - training data
            addParameter(p,'NumObservationsToSample',0);
            addParameter(p,'QueryPoints',[]);
            addParameter(p,'UseParallel',false);

            % Allowing additional Name-Value pairs specific to custom models
            % here so that a more specific error message can be thrown when parsing.
            addParameter(p,'CategoricalPredictors',[]);
            addParameter(p,'OutputColumns',[]);

            parse(p,model,features,labels,varargin{:});
            X = p.Results.Data;
            classif = true;

            % If third argument is a char, its a parameter name else it is Data
            % Data had already been parsed, so remove it from the varargin
            % list
            if(nargin>3 && ~ischar(varargin{1}))
                varargin = varargin(2:end);
            end
            [pd,x,y] = responsible.impl.partialDependence(model,features,...
                labels,X,classif,varargin{:});
            varargout{1} = x;
            varargout{2} = y;
        end

        function [AX] = plotPartialDependence(model,features,labels,varargin)
            %PLOTPARTIALDEPENDENCE Partial Dependence Plot for 1-D or 2-D visualization
            %   plotPartialDependence(MODEL,VARS,LABELS) plots the partial dependence
            %   PD between the predictor variables listed in VARS and the scores for
            %   the classes specified in LABELS by using the MultinomialRegression
            %   model MODEL and by averaging over the data used in fitting the model.
            %   You can specify one or two variables using VARS. VARS is either column
            %   indices of predictors or names of predictors. The names must match the
            %   entries in MODEL.PredictorNames. You can specify one or multiple classes
            %   using LABELS. The values in LABELS must match and have the same data type
            %   as the class names in MODEL.ClassNames. LABELS must be a single label if
            %   VARS specifies two predictors.
            %
            %   plotPartialDependence(...,DATA) specifies the data to be used for
            %   averaging. DATA is a matrix or table of data to be used in place
            %   of the data used in fitting the model.
            %
            %   plotPartialDependence(...,Name,Value) specifies additional options using
            %   one or more name-value pair arguments:
            %
            %   'Conditional'               'none' (default) instructs the method to create
            %                               a partial dependence plot with no conditioning,
            %                               'absolute' to create an individual conditional
            %                               expectation (ICE) plot, and 'centered' to
            %                               create an ICE plot with centered data.
            %
            %   'NumObservationsToSample'   An integer specifying the number
            %                               of rows to sample at random from the dataset
            %                               (either the data specified by the argument DATA
            %                               or the data used to train MODEL). The default
            %                               is to use all rows.
            %
            %   'QueryPoints'               The points at which to calculate the
            %                               partial dependence. With one predictor, the
            %                               value of 'QueryPoints' must be a column
            %                               vector. With two predictors, the value of
            %                               'QueryPoints' must be a 1x2 cell array
            %                               containing a separate vector for each
            %                               predictor. The default is to use 100 points
            %                               equally spaced across the range of the
            %                               predictor.
            %
            %   'UseParallel'               Can be true to instruct the method to perform
            %                               the averaging calculations in parallel (using
            %                               parfor), or false (default) to instruct the
            %                               method not to use parallel computations.
            %
            % 	'Parent'                    Instructs the method to create a partial
            %                               dependence plot using the axes with handle
            %                               specified by the value of this parameter.
            %
            %   AX = plotPartialDependence(...) returns a handle AX to the axes of the
            %   plot.
            %
            %   Examples:
            %       load census1994
            %       X = adultdata(:,{'age','workClass','education_num','marital_status','race',...
            %           'sex','capital_gain','capital_loss','hours_per_week','salary'});
            %       f = fitmnr(X,'salary');
            %
            %       plotPartialDependence(f,'age','<=50K');
            %       plotPartialDependence(f,{'age','workClass'},'>50K');
            %       plotPartialDependence(f,'age', {'<=50K' '>50K'})
            %
            %       % Obtain optional output Axes handle
            %       ax = plotPartialDependence(f,'age','<=50K');
            %
            %       % With additional Data
            %       plotPartialDependence(f,'age','<=50K',X);
            %       plotPartialDependence(f,1,'<=50K',X);
            %
            %       % With optional name-value pairs
            %       plotPartialDependence(f,1,'<=50K',X,'NumObservationsToSample',100);
            %       plotPartialDependence(f,1,'<=50K',X,'UseParallel',true);
            %       plotPartialDependence(f,1,'<=50K',X,'UseParallel',true,'Conditional','none');
            %
            %       % Plot the Individual Conditional Expectation
            %       plotPartialDependence(f,1,'<=50K',X,'Conditional','absolute');
            %
            %       % Provide alternative query points
            %       xi = linspace(min(X{:,'age'}),min(X{:,'age'}))';
            %       plotPartialDependence(f,1,'<=50K','QueryPoints',xi);
            %
            %       % Provide alternative query points
            %       xi = linspace(min(X{:,'age'}),max(X{:,'age'}))';
            %       plotPartialDependence(f,1,'<=50K','QueryPoints',xi);
            %
            %       xi = cell(1,2);
            %       xi{1} = linspace(min(X{:,1}),max(X{:,1}))';
            %       xi{2} = linspace(min(X{:,9}),max(X{:,9}))';
            %       plotPartialDependence(f,[1,9],'<=50K','QueryPoints',xi);
            %
            %   See also plotSlice, partialDependence, predict
            features = convertStringsToChars(features);
            [varargin{:}] = convertStringsToChars(varargin{:});
            p = inputParser;
            addRequired(p,'Model');
            addRequired(p,'Var');
            addRequired(p,'Label');
            addOptional(p,'Data',getData(model)); % Default - training data
            addParameter(p,'Conditional',{'none' 'absolute' 'centered'});
            addParameter(p,'NumObservationsToSample',0);
            addParameter(p,'QueryPoints',[]);
            addParameter(p,'UseParallel',false);
            addParameter(p,'ParentAxisHandle',[]);

            % Allowing additional Name-Value pairs specific to custom models
            % here so that a more specific error message can be thrown when parsing.
            addParameter(p,'CategoricalPredictors',[]);
            addParameter(p,'OutputColumns',[]);

            parse(p,model,features,labels,varargin{:});
            X = p.Results.Data;
            classif = true;

            % If third argument is a char, its a parameter name else it is Data
            % Data has already been parsed, so remove it from the varargin
            % list
            if(nargin>3 && ~ischar(varargin{1}))
                % Pass everything but the first argument(Data)to compact method
                varargin = varargin(2:end);
            end
            ax = responsible.impl.plotPartialDependence(model,features,...
                labels,X,classif,varargin{:});
            if(nargout > 0)
                AX = ax;
            end
        end

        function hOut = plotSlice(varargin)
            %PLOTSLICE Plot slices through the fitted regression surface
            %   PLOTSLICE(MNR) creates a new figure containing a series of
            %   plots, each representing a slice through the regression surface
            %   predicted by the MultinomialRegression MNR for
            %   MNR.ClassNames(1). For each plot, the surface slice is shown
            %   as a function of a single predictor variable, with the other
            %   predictor variables held constant, for a single class.
            %
            %   PLOTSLICE(MNR, XNEW) creates a a new figure using the
            %   predictor data in XNEW. XNEW must be a table or a
            %   double/single matrix, and must be able to be passed to the
            %   MNR PREDICT method. The predictor data in XNEW are used to 
            %   set the axes limits when generating the regression surface.
            %
            %   PLOTSLICE(MNR, PLOTTYPE) creates a plot through the
            %   regression surface in a style determined by PLOTTYPE.
            %   PLOTTYPE can be 'interactive' (default), 'histogram',
            %   or 'stackedhist'. 'interactive' generates the figure
            %   described above. 'histogram' generates a plot of
            %   histograms through the regression surface for each class
            %   overlaid on eachother. 'stackedhist' generates a plot of
            %   histograms through the regression surface for each class
            %   stacked on top of eachother. For both 'histogram' and 'stackedhist',
            %   the surfaces are generated by varying MNR.PredictorNames(1)
            %   and holding the remainder of the predictors fixed.
            %
            %   PLOTSLICE(MNR, PLOTTYPE, XNEW) creates one of the plots above
            %   in the style of PLOTTYPE using the predictor data in XNEW. 
            %   XNEW must be a table or a double/single matrix, and must 
            %   be able to be passed to the MNR PREDICT method. The predictor 
            %   data in XNEW are used to set the axes limits when generating 
            %   the regression surface.
            %
            %   PLOTSLICE(__, Name1=Val1, ...) specifies optional
            %   name-value arguments:
            %       'ClassToPlot' - Which class among MNR.ClassNames to show
            %                       the regression surface for. This is only
            %                       relevant when PLOTTYPE is 'interactive'.
            %                       By default, this is MNR.ClassNames(1).
            %
            %       'PredictorToVary' - Which predictor among MNR.PredictorNames
            %                           to vary when creating the surface.
            %                           This is only relevant when PLOTTYPE
            %                           is 'histogram' or 'stackedhist'. By
            %                           default, this is MNR.PredictorNames(1).
            %
            %   The PLOTTYPE argument can be followed by name-value arguments
            %   to specify additional properties of the objects in the plot.
            %   If PLOTTYPE is 'histogram', you can specify Histogram object
            %   properties. Otherwise, if PLOTTYPE is 'stackedhist', you can
            %   specify Bar object properties.
            %
            %   PLOTSLICE(FIG, __) plots into the figure FIG. FIG must be 
            %   created by the uifigure function, and is only supported with
            %   PLOTTYPE is 'interactive'.
            %
            %   PLOTSLICE(AX,__) plots into the specified axes AX. PLOTTYPE
            %   must be 'histogram' or 'stackedhist'.
            %
            %   H = PLOTSLICE(__) returns a handle H to the Figure (if PLOTTYPE
            %   is 'interactive'), Histogram objects (if PLOTTYPE is 'histogram'),
            %   or Bar objects (if PLOTTYPE is 'stackedhist') in the plot.
            %
            %
            %   Example: Plot an interactive figure through the regression
            %            surface
            %   load fisheriris
            %   mnr = fitmnr(meas, species);
            %   plotSlice(mnr);
            %
            %   See also plotResiduals, plotPartialDependence.
            [parent, model, varargin] = MultinomialRegression.parseParentAndModel(varargin{:});
            types = ["interactive", "histogram", "stackedhist"];
            numVararg = numel(varargin);
            if mod(numVararg, 2) == 0
                % Even - matched NV args, or plottype, predictor data,
                % and matched NV args.
                if numVararg > 0
                    if internal.stats.isString(varargin{1}) && ismember(varargin{1}, types)
                        % plottype, predictor data, and NV args
                        plottype = varargin{1};
                        preddata = varargin{2};
                        varargin = varargin(3:end);
                    else
                        % NV args 
                        plottype = 'interactive';
                        preddata = [];
                    end
                else
                    plottype = 'interactive';
                    preddata = [];
                end
            else
                % Odd - either plottype and NV args or predictor data and
                % NV args
                if internal.stats.isString(varargin{1})
                    plottype = varargin{1};
                    plottype = validatestring(plottype, types, 'plotSlice', 'plottype');
                    preddata = [];
                else
                    % Let the predict method validate the predictor data
                    % later on
                    plottype = 'interactive';
                    preddata = varargin{1};
                end
                varargin = varargin(2:end);
            end
            if strcmp(plottype, 'interactive')
                h = model.plotInteractive(model, preddata, parent, varargin{:});
            else
                h = model.plotHistogram(model, parent, plottype, preddata, varargin{:});
            end
            if nargout == 1
                hOut = h;
            end
        end

        function hOut = plotResiduals(varargin)
            %PLOTRESIDUALS Plot residuals of fitted model
            %   PLOTRESIDUALS(MNR) plots the deviance residuals from
            %   MultinomialRegression MNR in a histogram plot.
            %
            %   PLOTRESIDUALS(MNR, PLOTTYPE) plots the deviance residuals
            %   from MultinomialRegression MNR in a plot of type PLOTTYPE.
            %   Valid values for PLOTTYPE are:
            %
            %      'caseorder'     residuals vs. case (row) order
            %      'histogram'     histogram (default)
            %      'lagged'        residual vs. lagged residual (r(t) vs. r(t-1))
            %      'probability'   normal probability plot
            %      'symmetry'      symmetry plot
            %
            %
            %   plotResiduals(__, NAME1=VAL1, ...) specifies optional
            %   name-value arguments:
            %      'ResidualType' - The type of residual to show in the
            %                        plot. The supported types are:
            %                            'raw'      - Observed minus fitted values
            %                            'Pearson'  - Residuals divided by the estimated
            %                                         standard deviation of the response.
            %                            'deviance' - Residuals based on the
            %                                         contribution of each
            %                                         observation to the deviance.
            %                                         (default)
            %
            %      'ClassToPlot' - Which class among MNR.ClassNames to
            %                      show the plot of residuals for. This
            %                      is only relevant when 'ResidualType'
            %                      is not 'deviance'. The raw and Pearson
            %                      residuals are computed per-class, rather
            %                      than per the whole model for deviance
            %                      residuals. By default, MNR.ClassNames(1)
            %                      is chosen.
            %
            %   The PLOTTYPE argument can be followed by name-value arguments
            %   to specify additional properties of the objects in the plot.
            %   If PLOTTYPE is 'histogram', you can specify Patch object
            %   properties. Otherwise, if PLOTTYPE is not 'probability',
            %   you can specify Line object properties.
            %
            %   PLOTRESIDUALS(AX,__) plots into the specified axes AX.
            %
            %   H = PLOTRESIDUALS(__) returns a handle to the patch (if
            %   PLOTYPE is 'histogram') or lines (for other PLOTTYPES) in
            %   the plot.
            %
            %   Example: Plot the residuals against the case order
            %   load fisheriris
            %   mnr = fitmnr(meas, species);
            %   plotResiduals(mnr, 'caseorder');
            %
            %   See also plotSlice, plotPartialDependence, testDeviance
            [ax, model, varargin] = MultinomialRegression.parseParentAndModel(varargin{:});

            if mod(numel(varargin), 2) == 0
                % Even - matched NV args
                % Use default plot type 'histogram'
                plottype = 'histogram';
            else
                % Odd - plottype + NV args
                plottype = varargin{1};
                varargin = varargin(2:end);
            end

            args = {'ClassToPlot', 'ResidualType'};
            defs = {model.ClassNames(1), 'deviance'};
            [className, residType, wasSet, varargin] = internal.stats.parseArgs(args, defs, varargin{:});
            model.validateClassName(className, 'scalar', 'ClassToPlot');
            residType = validatestring(residType, {'raw','pearson','deviance'}, 'plotResiduals', 'ResidualType');
            plottype = validatestring(plottype, {'histogram', 'caseorder', 'lagged', 'probability', 'symmetry'}, ...
                'plotResiduals', 'plottype');
            className = classreg.learning.internal.ClassLabel(className);
            [tf, pos] = ismember(className, model.PrivClassNames);
            pos = pos(tf);
            if ~tf
                error(message('stats:MultinomialRegression:UnrecognizedClass', char(className)));
            end

            if isempty(ax)
                ax = gca();
            end
            if strcmpi(residType, 'deviance')
                % Deviance residuals are scalar
                pos = 1;
                if wasSet.ClassToPlot
                    % Ignored for deviance residuals, as they are for the
                    % whole model, unlike other residuals which are per
                    % class
                    warning(message('stats:MultinomialRegression:ClassToPlotIgnored'))
                end
            end
            h = classreg.regr.modelutils.plotResiduals(model, plottype, varargin{:}, ...
                'ResidualType', residType, 'ResidualIndex', pos, 'Parent', ax);
            if nargout > 0
                hOut = h;
            end
        end
    end

    methods (Hidden)
        function disp(model)
            %DISP Display a MultinomialRegression.
            %   DISP(MODEL) displays the MultinomialRegression MODEL.
            %
            %   See also MultinomialRegression.

            % Adjust spacing based on Command Window spacing setting
            spacing = settings().matlab.commandwindow.DisplayLineSpacing;
            dispStr = '';
            if strcmp(spacing.ActiveValue, 'loose')
                dispStr = '\n';
            end

            % Get header line
            dispStr = [dispStr, getString(message('stats:MultinomialRegression:display_MultinomialRegression', model.ModelType))];

            % Get coefficient table
            dispStr = [dispStr, newline, newline, evalc('disp(model.Coefficients)')];

            % Get DFE, number of obs
            dispStr = [dispStr, getString(message('stats:GeneralizedLinearModel:display_NumObsDFE', model.NumObservations, model.DFE))];
            [p, stat] = testDeviance(model);
            p = num2str(p);
            stat = num2str(stat);

            % Get the Chi^2/F test results
            if model.EstimateDispersion
                disperStr = getString(message('stats:GeneralizedLinearModel:display_EstimatedDispersion', num2str(model.Dispersion)));
                testStr = getString(message('stats:GeneralizedLinearModel:display_FTest', stat, p));
            else
                disperStr = getString(message('stats:GeneralizedLinearModel:display_Dispersion', num2str(model.Dispersion)));
                testStr = getString(message('stats:GeneralizedLinearModel:display_ChiTest', stat, p));
            end
            dispStr = [dispStr, disperStr, testStr];
            fprintf(dispStr);
        end

        %%% OTHER INHERITED METHODS FROM PARAMETRICREGRESSION %%%
        function t = title(model)
            strLHS = classreg.regr.CompactGeneralizedLinearModel.linkstr(model.Formula.Link,model.ResponseName);
            strFunArgs = internal.stats.strCollapse(model.Formula.PredictorNames,',');
            t = sprintf('%s = mnr(%s)',strLHS,strFunArgs);
        end

        function out = varianceParam(model)
            out = model.Dispersion;
        end
    end
end