www.gusucode.com > stats 源码程序 matlab案例代码 > stats/PlottheCountEstimatesExample.m

    %% Plot the Count Estimates  

% Copyright 2015 The MathWorks, Inc.


%% 
% Create sample data with one predictor variable and a categorical response
% variable with three categories. 
x = [-3 -2 -1 0 1 2 3]';
Y = [1 11 13; 2 9 14; 6 14 5; 5 10 10;...
		 5 14 6; 7 13 5; 8 11 6];
[Y x] 

%%
% There are observations on seven different values of the predictor
% variable |x| . The response variable |Y| has three categories and the data
% shows how many of the 25 individuals are in each category of |Y| for each
% observation of |x|. For example, when |x| is -3, 1 of 25 individuals is
% observed in category 1, 11 observed in category 2, and 13 observed in
% category 3. Similarly, when |x| is 1, 5 of the individuals are observed
% in category 1, 14 are observed in category 2, and 6 are observed in
% category 3.

%% 
% Plot the number in each category versus the |x| values, on a stacked bar
% graph. 
bar(x,Y,'stacked'); 
ylim([0 25]);     

%% 
% Fit a nominal model for the individual response category probabilities,
% with separate slopes on the single predictor variable, |x|, for each
% category.
betaHatNom = mnrfit(x,Y,'model','nominal',...
    'interactions','on') 

%%
% The first row of |betaHatOrd| contains the intercept terms for the first
% two response categories. The second row contains the slopes. |mnrfit|
% accepts the third category as the reference category and hence assumes
% the coefficients for the third category are zero.  

%% 
% Compute the predicted probabilities for the three response categories. 
xx = linspace(-4,4)';
piHatNom = mnrval(betaHatNom,xx,'model','nominal',...
    'interactions','on'); 

%%
% The probability of being in the third category is simply 1 - P($y$ = 1)
% - P($y$ = 2).  

%% 
% Plot the estimated cumulative number in each category on the bar graph. 
line(xx,cumsum(25*piHatNom,2),'LineWidth',2);    

%%
% The cumulative probability for the third category is always 1.  

%% 
% Now, fit a "parallel" ordinal model for the cumulative response category
% probabilities, with a common slope on the single predictor variable, |x|,
% across all categories: 
betaHatOrd = mnrfit(x,Y,'model','ordinal',...
    'interactions','off') 

%%
% The first two elements of |betaHatOrd| are the intercept terms for the
% first two response categories. The last element of |betaHatOrd| is the
% common slope.  

%% 
% Compute the predicted cumulative probabilities for the first two response
% categories. The cumulative probability for the third category is always
% 1.
piHatOrd = mnrval(betaHatOrd,xx,'type','cumulative',...
    'model','ordinal','interactions','off');  

%% 
% Plot the estimated cumulative number on the bar graph of the observed
% cumulative number. 
figure()
bar(x,cumsum(Y,2),'grouped'); 
ylim([0 25]);
line(xx,25*piHatOrd,'LineWidth',2);