- 在线时间
- 326 小时
- 最后登录
- 2019-9-18
- 注册时间
- 2014-8-5
- 听众数
- 36
- 收听数
- 9
- 能力
- 0 分
- 体力
- 4485 点
- 威望
- 0 点
- 阅读权限
- 60
- 积分
- 1854
- 相册
- 1
- 日志
- 0
- 记录
- 0
- 帖子
- 996
- 主题
- 413
- 精华
- 0
- 分享
- 3
- 好友
- 98
升级   85.4% TA的每日心情 | 开心 2019-9-18 21:55 |
---|
签到天数: 258 天 [LV.8]以坛为家I
群组: 2015国赛冲刺 群组: 2016美赛公益课程 群组: 国赛讨论 群组: 第三届数模基础实训 群组: Matlab讨论组 |
- function [p,anovatab,stats] = anova1(x,group,displayopt,extra)
- %ANOVA1 One-way analysis of variance (ANOVA).
- % ANOVA1 performs a one-way ANOVA for comparing the means of two or more
- % groups of data. It returns the p-value for the null hypothesis that the
- % means of the groups are equal.
- %
- % P = ANOVA1(M) for a matrix M treats each column as a separate group,
- % and determines whether the population means of the columns are equal.
- % This form of ANOVA1 is appropriate when each group has the same number
- % of elements (balanced ANOVA).
- %
- % P = ANOVA1(V,GROUP) groups elements in the vector V according to values
- % in the grouping variable GROUP. GROUP must be a categorical variable,
- % numeric vector, logical vector, string array, or cell array of strings
- % with one group name for each element of X. X values corresponding to
- % the same value of GROUP are placed in the same group.
- %
- % P = ANOVA1(M,GROUP) accepts a character array or cell array of strings,
- % with one group name for each column of M. Columns with the same group
- % name are treated as part of the same group.
- %
- % P = ANOVA1(X,GROUP,DISPLAYOPT) controls the display. DISPLAYOPT can be
- % 'on' (the default) to display figures containing a standard one-way
- % anova table and a boxplot, or 'off' to omit these displays. Note that
- % the notches in the boxplot provide a test of group medians (see HELP
- % BOXPLOT), and this is not the same as the F test for different means in
- % the anova table. X can be either a vector or matrix. If X is a matrix
- % and there are no group names, specify GROUP as [].
- %
- % [P,ANOVATAB] = ANOVA1(...) returns the ANOVA table values as the
- % cell array ANOVATAB.
- %
- % [P,ANOVATAB,STATS] = ANOVA1(...) returns an additional structure
- % of statistics useful for performing a multiple comparison of means
- % with the MULTCOMPARE function.
- %
- % See also ANOVA2, ANOVAN, BOXPLOT, MANOVA1, MULTCOMPARE.
- % Reference: Robert V. Hogg, and Johannes Ledolter, Engineering Statistics
- % Macmillan 1987 pp. 205-206.
- % Copyright 1993-2013 The MathWorks, Inc.
- narginchk(1,4);
- classical = 1;
- nargs = nargin;
- if (nargin>0 && strcmp(x,'kruskalwallis'))
- % Called via kruskalwallis function, adjust inputs
- classical = 0;
- if (nargin >= 2), x = group; group = []; end
- if (nargin >= 3), group = displayopt; displayopt = []; end
- if (nargin >= 4), displayopt = extra; end
- nargs = nargs-1;
- end
- if (nargs < 2), group = []; end
- if (nargs < 3), displayopt = 'on'; end
- % Note: for backwards compatibility, accept 'nodisplay' for 'off'
- willdisplay = ~(strcmp(displayopt,'nodisplay') | strcmp(displayopt,'n') ...
- | strcmp(displayopt,'off'));
- % Convert group to cell array from character array, make it a column
- if (ischar(group) && ~isempty(group))
- group = cellstr(group);
- end
- if (size(group, 1) == 1)
- group = group';
- end
- % If the input is a matrix, it may not be balanced if it contains NaNs or
- % if there are repeated grouping values, so turn it into a vector.
- needvector = false;
- if ~isvector(x)
- if (any(isnan(x(:))))
- needvector = true;
- elseif ~isempty(group) && length(group)~=length(unique(group))
- needvector = true;
- end
- end
- % If X is a matrix with NaNs, convert to vector form.
- if ~isvector(x)
- if needvector
- [n,m] = size(x);
- x = x(:);
- gi = reshape(repmat((1:m), n, 1), n*m, 1);
- if isempty(group) % no group names
- group = gi;
- elseif (size(group,1) == m)
- group = group(gi,:);
- else
- error(message('stats:anova1:InputSizeMismatch'));
- end
- end
- elseif ~isempty(group) && (size(group,1) ~= length(x))
- error(message('stats:anova1:InputSizeMismatch'));
- end
- % If X is a matrix GROUP is provided with the correct size, use GROUP
- % to define groups and to label boxes
- if (~isempty(group) && (length(x) < numel(x)) ...
- && (size(x,2) == size(group,1)))
- named = 1;
- [gid,gnames] = grp2idx(group);
- gnames = gnames(gid);
- grouped = 0;
- else
- named = 0;
- gnames = [];
- grouped = ~isempty(group);
- end
- if (grouped)
- % Single data vector and a separate grouping variable
- x = x(:);
- lx = length(x);
- if (lx ~= numel(x))
- error(message('stats:anova1:VectorRequired'))
- end
- nonan = ~isnan(x);
- x = x(nonan);
- % Convert group to indices 1,...,g and separate names
- group = group(nonan,:);
- [groupnum, gnames] = grp2idx(group);
- named = 1;
- % Remove NaN values
- nonan = ~isnan(groupnum);
- if (~all(nonan))
- groupnum = groupnum(nonan);
- x = x(nonan);
- end
- lx = length(x);
- xorig = x; % use uncentered version to make M
- groupnum = groupnum(:);
- maxi = size(gnames, 1);
- if isa(x,'single')
- xm = zeros(1,maxi,'single');
- else
- xm = zeros(1,maxi);
- end
- countx = xm;
- if (classical)
- mu = mean(x);
- x = x - mu; % center to improve accuracy
- xr = x;
- else
- [xr,tieadj] = tiedrank(x);
- end
-
- for j = 1:maxi
- % Get group sizes and means
- k = find(groupnum == j);
- lk = length(k);
- countx(j) = lk;
- xm(j) = mean(xr(k)); % column means
- end
- gm = mean(xr); % grand mean
- df1 = sum(countx>0) - 1; % Column degrees of freedom
- df2 = lx - df1 - 1; % Error degrees of freedom
- xc = xm - gm; % centered
- xc(countx==0) = 0;
- RSS = dot(countx, xc.^2); % Regression Sum of Squares
- else
- % Data in matrix form, no separate grouping variable
- [r,c] = size(x);
- lx = r * c;
- if (classical)
- xr = x;
- mu = mean(xr(:));
- xr = xr - mu; % center to improve accuracy
- else
- [xr,tieadj] = tiedrank(x(:));
- xr = reshape(xr, size(x));
- end
- countx = repmat(r, 1, c);
- xorig = x; % save uncentered version for boxplot
- xm = mean(xr); % column means
- gm = mean(xm); % grand mean
- df1 = c-1; % Column degrees of freedom
- df2 = c*(r-1); % Error degrees of freedom
- RSS = r*(xm - gm)*(xm-gm)'; % Regression Sum of Squares
- end
- TSS = (xr(:) - gm)'*(xr(:) - gm); % Total Sum of Squares
- SSE = TSS - RSS; % Error Sum of Squares
- if (df2 > 0)
- mse = SSE/df2;
- else
- mse = NaN;
- end
- if (classical)
- if (SSE~=0)
- F = (RSS/df1) / mse;
- p = fpval(F,df1,df2); % Probability of F given equal means.
- elseif (RSS==0) % Constant Matrix case.
- F = NaN;
- p = NaN;
- else % Perfect fit case.
- F = Inf;
- p = 0;
- end
- else
- F = (12 * RSS) / (lx * (lx+1));
- if (tieadj > 0)
- F = F / (1 - 2 * tieadj/(lx^3-lx));
- end
- p = chi2pval(F,df1);
- end
- Table=zeros(3,5); %Formatting for ANOVA Table printout
- Table(:,1)=[ RSS SSE TSS]';
- Table(:,2)=[df1 df2 df1+df2]';
- Table(:,3)=[ RSS/df1 mse Inf ]';
- Table(:,4)=[ F Inf Inf ]';
- Table(:,5)=[ p Inf Inf ]';
- colheads = {getString(message('stats:anova1:ColHeadSource')), getString(message('stats:anova1:ColHeadSS')), getString(message('stats:anova1:ColHeadDf')), getString(message('stats:anova1:ColHeadMS')), getString(message('stats:anova1:ColHeadF')), getString(message('stats:anova1:ColHeadProbGtF'))};
- if (~classical)
- colheads{5} = getString(message('stats:anova1:ColHeadChisq'));
- colheads{6} = getString(message('stats:anova1:ColHeadProbGtChisq'));
- end
- rowheads = {getString(message('stats:anova1:RowHeadColumns')), getString(message('stats:anova1:RowHeadError')), getString(message('stats:anova1:RowHeadTotal'))};
- if (grouped)
- rowheads{1} = getString(message('stats:anova1:RowHeadGroups'));
- end
- % Create cell array version of table
- atab = num2cell(Table);
- for i=1:size(atab,1)
- for j=1:size(atab,2)
- if (isinf(atab{i,j}))
- atab{i,j} = [];
- end
- end
- end
- atab = [rowheads' atab];
- atab = [colheads; atab];
- if (nargout > 1)
- anovatab = atab;
- end
- % Create output stats structure if requested, used by MULTCOMPARE
- if (nargout > 2)
- if ~isempty(gnames)
- stats.gnames = gnames;
- else
- stats.gnames = strjust(num2str((1:length(xm))'),'left');
- end
- stats.n = countx;
- if (classical)
- stats.source = 'anova1';
- stats.means = xm + mu;
- stats.df = df2;
- stats.s = sqrt(mse);
- else
- stats.source = 'kruskalwallis';
- stats.meanranks = xm;
- stats.sumt = 2 * tieadj;
- end
- end
- if (~willdisplay), return; end
- digits = [-1 -1 0 -1 2 4];
- if (classical)
- wtitle = getString(message('stats:anova1:OnewayANOVA'));
- ttitle = getString(message('stats:anova1:ANOVATable'));
- else
- wtitle = getString(message('stats:anova1:KruskalWallisOnewayANOVA'));
- ttitle = getString(message('stats:anova1:KruskalWallisANOVATable'));
- end
- tblfig = statdisptable(atab, wtitle, ttitle, '', digits);
- set(tblfig,'tag','table');
- f1 = figure('pos',get(gcf,'pos') + [0,-200,0,0],'tag','boxplot');
- ax = axes('Parent',f1);
- if (~grouped)
- boxplot(ax,xorig,'notch','on');
- else
- boxplot(ax,xorig,groupnum,'notch','on');
- h = get(ax,'XLabel');
- set(h,'String',getString(message('stats:anova1:GroupNumber')));
- end
- % If there are group names, use them
- if ~isempty(gnames)
- h = get(ax,'XLabel');
- if (named)
- set(h,'String','');
- end
- set(ax, 'xtick', (1:size(gnames,1)), 'xticklabel', gnames);
- end
复制代码
|
|