判别分析与主成分分析

实验目的

　　（1）掌握判别分析、主成分分析。

　　（2）会用判别分析、主成分分析对实际问题进行分析。

实验要求

 　　实验步骤要有模型建立，模型求解、结果分析。

实验内容

　　（1）银行的贷款部门需要判别每个客户的信用好坏（是否未履行还贷责任），以决定是否给予贷款。可以根据贷款申请人的年龄（X1）、受教育程度（X2）、现在所从事工作的年数（X3）、未变更住址的年数（X4）、收入（X5）、负债收入比例（X6）、信用卡债务（X7）、其它债务（X8）等来判断其信用情况。下表是从某银行的客户资料中抽取的部分数据，和某客户的如上情况资料为（53，1，9，18，50，11.20，2.02，3.58），根据样本资料分别用马氏距离判别法、线性判别法、二次判别法对其进行信用好坏的判别。 

目前信用好坏	客户序号	X1	X2	X3	X4	X5	X6	X7	X8
已履行还贷责任	1	23	1	7	2	31	6.60	0.34	1.71
	2	34	1	17	3	59	8.00	1.81	2.91
	3	42	2	7	23	41	4.60	0.94	.94
	4	39	1	19	5	48	13.10	1.93	4.36
	5	35	1	9	1	34	5.00	0.40	1.30
未履行还贷责任	6	37	1	1	3	24	15.10	1.80	1.82
	7	29	1	13	1	42	7.40	1.46	1.65
	8	32	2	11	6	75	23.30	7.76	9.72
	9	28	2	2	3	23	6.40	0.19	1.29
	10	26	1	4	3	27	10.50	2.47	.36

　　（2）在某中学随机抽取某年级30名学生，测量其身高（X1）、体重（X2）、胸围（X3）和坐高（X4），数据如表（30名中学生身体四项指标数据），试对这30名中学生身体四项指标数据做主成分分析。

 实验步骤

 

组统计
group		平均值	标准 偏差	有效个案数（成列）
				未加权	加权
1.00	X1	34.6000	7.23187	5	5.000
	X2	1.2000	.44721	5	5.000
	X3	11.8000	5.76194	5	5.000
	X4	6.8000	9.17606	5	5.000
	X5	42.6000	11.28273	5	5.000
	X6	7.4600	3.43045	5	5.000
	X7	1.0840	.75580	5	5.000
	X8	2.2440	1.39622	5	5.000
2.00	X1	30.4000	4.27785	5	5.000
	X2	1.4000	.54772	5	5.000
	X3	6.2000	5.44977	5	5.000
	X4	3.2000	1.78885	5	5.000
	X5	38.2000	21.94767	5	5.000
	X6	12.5400	6.90312	5	5.000
	X7	2.7360	2.92821	5	5.000
	X8	2.9680	3.81647	5	5.000
总计	X1	32.5000	6.02310	10	10.000
	X2	1.3000	.48305	10	10.000
	X3	9.0000	6.05530	10	10.000
	X4	5.0000	6.51494	10	10.000
	X5	40.4000	16.61459	10	10.000
	X6	10.0000	5.79463	10	10.000
	X7	1.9100	2.19609	10	10.000
	X8	2.6060	2.73598	10	10.000

 　　由上表得到下列式子：

　　　　　　　　　　　　　　 　　　　　　　　　　

协方差矩阵
group		X1	X2	X3	X4	X5	X6	X7	X8
1.00 S1	X1	52.300	1.850	11.900	41.900	33.300	3.080	2.645	1.269
	X2	1.850	.200	-1.200	4.050	-.400	-.715	-.036	-.326
	X3	11.900	-1.200	33.200	-17.800	52.900	17.040	4.011	7.541
	X4	41.900	4.050	-17.800	84.200	1.900	-10.035	.231	-4.857
	X5	33.300	-.400	52.900	1.900	127.300	18.755	7.805	9.687
	X6	3.080	-.715	17.040	-10.035	18.755	11.768	1.974	4.701
	X7	2.645	-.036	4.011	.231	7.805	1.974	.571	.876
	X8	1.269	-.326	7.541	-4.857	9.687	4.701	.876	1.949
2.00 S2	X1	18.300	-.200	-4.100	1.900	11.400	16.255	2.732	5.144
	X2	-.200	.300	.150	.650	5.400	1.155	.620	1.269
	X3	-4.100	.150	29.700	.200	91.200	8.415	7.896	10.551
	X4	1.900	.650	.200	3.200	25.700	10.640	4.406	5.723
	X5	11.400	5.400	91.200	25.700	481.700	114.065	58.751	78.620
	X6	16.255	1.155	8.415	10.640	114.065	47.653	18.599	23.028
	X7	2.732	.620	7.896	4.406	58.751	18.599	8.574	10.411
	X8	5.144	1.269	10.551	5.723	78.620	23.028	10.411	14.565

　　S₁与S₂见上表。

　　样本协方差矩阵由上表得到。

汇聚组内矩阵a
		X1	X2	X3	X4	X5	X6	X7	X8
协方差	X1	35.300	.825	3.900	21.900	22.350	9.668	2.688	3.206
	X2	.825	.250	-.525	2.350	2.500	.220	.292	.471
	X3	3.900	-.525	31.450	-8.800	72.050	12.728	5.953	9.046
	X4	21.900	2.350	-8.800	43.700	13.800	.303	2.319	.433
	X5	22.350	2.500	72.050	13.800	304.500	66.410	33.278	44.154
	X6	9.668	.220	12.728	.303	66.410	29.711	10.287	13.864
	X7	2.688	.292	5.953	2.319	33.278	10.287	4.573	5.644
	X8	3.206	.471	9.046	.433	44.154	13.864	5.644	8.257
相关性	X1	1.000	.278	.117	.558	.216	.299	.212	.188
	X2	.278	1.000	-.187	.711	.287	.081	.273	.328
	X3	.117	-.187	1.000	-.237	.736	.416	.496	.561
	X4	.558	.711	-.237	1.000	.120	.008	.164	.023
	X5	.216	.287	.736	.120	1.000	.698	.892	.881
	X6	.299	.081	.416	.008	.698	1.000	.883	.885
	X7	.212	.273	.496	.164	.892	.883	1.000	.918
	X8	.188	.328	.561	.023	.881	.885	.918	1.000
a. 协方差矩阵的*度为 8。

 

分类函数系数
	group
	1.00	2.00
X1	.340	.184
X2	94.070	126.660
X3	1.033	1.874
X4	-4.943	-6.681
X5	2.969	3.086
X6	13.723	17.182
X7	-10.994	-7.133
X8	-37.504	-49.116
(常量)	-118.693	-171.296
费希尔线性判别函数

 　　第一类所对应的判别函数：

　　

　　第二类所对应的判别函数：

　　

 　　3、判别效果验证，直接验证10个个体的原始数据

分类结果a
		group	预测组成员信息		总计
			1.00	2.00
原始	计数	1.00	5	0	5
		2.00	0	5	5
		未分组个案	1	0	1
	%	1.00	100.0	.0	100.0
		2.00	.0	100.0	100.0
		未分组个案	100.0	.0	100.0
a. 正确地对 100.0% 个原始已分组个案进行了分类。

　　第1类别的正确率100%；第2类别的正确类别100%。

　　最终结果，该客户为第1类“已履行还贷责任”。

 　　代码：

DISCRIMINANT
  /GROUPS=group(1 2)
  /VARIABLES=X1 X2 X3 X4 X5 X6 X7 X8
  /ANALYSIS ALL
  /SAVE=CLASS SCORES PROBS
  /PRIORS EQUAL
  /STATISTICS=MEAN STDDEV UNIVF BOXM COEFF CORR COV GCOV TABLE
  /CLASSIFY=NONMISSING POOLED.

　　马氏距离判别法和二次判别法的详细步骤依次类推，下面使用MATLAB编程应用马氏距离判别法，线性判别法和二次判别法求解题1。

　　MATLAB代码：

% 样本数据
training=[23,1,7,2,31,6.60,0.34,1.71
34,1,17,3,59,8.00,1.81,2.91
42,2,7,23,41,4.60,0.94,.94
39,1,19,5,48,13.10,1.93,4.36
35,1,9,1,34,5.00,0.40,1.30
37,1,1,3,24,15.10,1.80,1.82
29,1,13,1,42,7.40,1.46,1.65
32,2,11,6,75,23.30,7.76,9.72
28,2,2,3,23,6.40,0.19,1.29
26,1,4,3,27,10.50,2.47,.36];
% 用于构造判别函数的样本数据矩阵
group=[1,1,1,1,1,2,2,2,2,2]\';
% 参数group是与training相应的分组变量
sample=[53.00,1.00,9.00,18.00,50.00,11.20,2.02,3.58];
% 使用线性判别法分类
[class2,err2]=classify(sample,training,group,\'linear\')
% 使用二次判别法分类
[class3,err3]=classify(sample,training,group,\'diagQuadratic\') 
% 使用马氏距离判别法分类
[class1,err1]=classify(sample,training,group,\'mahalanobis\')

　　运行示例：

>> exp64

class2 =

     1

err2 =

     0

class3 =

     1

err3 =

    0.1000

错误使用 classify (line 282)

The covariance matrix of each group in TRAINING must be positive definite.

　　在运行过程中，二次判别这里使用的对角的二次判别代替基于协方差的二次判别；显然线性判别与上一步的求解结果一致，而二次判别与线性判别到结果都是把样本分类为第1类“已履行还贷责任”一类。但是在调用MATLAB的classify时，本题的数据使用马氏距离判别出现了协方差矩阵非正定的错误。（这个错误的来源并非是协方差，而是在求协方差矩阵的逆的时候有一些数据出现了过小的情况）

　　本报告修改classify代码如下：

%             if any(s <= max(gsize(k),d) * eps(max(s)))
%                 error(message(\'stats:classify:BadQuadVar\'));
%             end

　　完整代码：

function [outclass, err, posterior, logp, coeffs] = classify(sample, training, group, type, prior)
%CLASSIFY Discriminant analysis.
%   CLASS = CLASSIFY(SAMPLE,TRAINING,GROUP) classifies each row of the data
%   in SAMPLE into one of the groups in TRAINING.  SAMPLE and TRAINING must
%   be matrices with the same number of columns.  GROUP is a grouping
%   variable for TRAINING.  Its unique values define groups, and each
%   element defines which group the corresponding row of TRAINING belongs
%   to.  GROUP can be a categorical variable, numeric vector, a string
%   array, or a cell array of strings.  TRAINING and GROUP must have the
%   same number of rows.  CLASSIFY treats NaNs or empty strings in GROUP as
%   missing values, and ignores the corresponding rows of TRAINING. CLASS
%   indicates which group each row of SAMPLE has been assigned to, and is
%   of the same type as GROUP.
%
%   CLASS = CLASSIFY(SAMPLE,TRAINING,GROUP,TYPE) allows you to specify the
%   type of discriminant function, one of \'linear\', \'quadratic\',
%   \'diagLinear\', \'diagQuadratic\', or \'mahalanobis\'.  Linear discrimination
%   fits a multivariate normal density to each group, with a pooled
%   estimate of covariance.  Quadratic discrimination fits MVN densities
%   with covariance estimates stratified by group.  Both methods use
%   likelihood ratios to assign observations to groups.  \'diagLinear\' and
%   \'diagQuadratic\' are similar to \'linear\' and \'quadratic\', but with
%   diagonal covariance matrix estimates.  These diagonal choices are
%   examples of naive Bayes classifiers.  Mahalanobis discrimination uses
%   Mahalanobis distances with stratified covariance estimates.  TYPE
%   defaults to \'linear\'.
%
%   CLASS = CLASSIFY(SAMPLE,TRAINING,GROUP,TYPE,PRIOR) allows you to
%   specify prior probabilities for the groups in one of three ways.  PRIOR
%   can be a numeric vector of the same length as the number of unique
%   values in GROUP (or the number of levels defined for GROUP, if GROUP is
%   categorical).  If GROUP is numeric or categorical, the order of PRIOR
%   must correspond to the ordered values in GROUP, or, if GROUP contains
%   strings, to the order of first occurrence of the values in GROUP. PRIOR
%   can also be a 1-by-1 structure with fields \'prob\', a numeric vector,
%   and \'group\', of the same type as GROUP, and containing unique values
%   indicating which groups the elements of \'prob\' correspond to. As a
%   structure, PRIOR may contain groups that do not appear in GROUP. This
%   can be useful if TRAINING is a subset of a larger training set.
%   CLASSIFY ignores any groups that appear in the structure but not in the
%   GROUPS array.  Finally, PRIOR can also be the string value \'empirical\',
%   indicating that the group prior probabilities should be estimated from
%   the group relative frequencies in TRAINING.  PRIOR defaults to a
%   numeric vector of equal probabilities, i.e., a uniform distribution.
%   PRIOR is not used for discrimination by Mahalanobis distance, except
%   for error rate calculation.
%
%   [CLASS,ERR] = CLASSIFY(...) returns ERR, an estimate of the
%   misclassification error rate that is based on the training data.
%   CLASSIFY returns the apparent error rate, i.e., the percentage of
%   observations in the TRAINING that are misclassified, weighted by the
%   prior probabilities for the groups.
%
%   [CLASS,ERR,POSTERIOR] = CLASSIFY(...) returns POSTERIOR, a matrix
%   containing estimates of the posterior probabilities that the j\'th
%   training group was the source of the i\'th sample observation, i.e.
%   Pr{group j | obs i}.  POSTERIOR is not computed for Mahalanobis
%   discrimination.
%
%   [CLASS,ERR,POSTERIOR,LOGP] = CLASSIFY(...) returns LOGP, a vector
%   containing estimates of the logs of the unconditional predictive
%   probability density of the sample observations, p(obs i) is the sum of
%   p(obs i | group j)*Pr{group j} taken over all groups.  LOGP is not
%   computed for Mahalanobis discrimination.
%
%   [CLASS,ERR,POSTERIOR,LOGP,COEF] = CLASSIFY(...) returns COEF, a
%   structure array containing coefficients describing the boundary between
%   the regions separating each pair of groups.  Each element COEF(I,J)
%   contains information for comparing group I to group J, defined using
%   the following fields:
%       \'type\'      type of discriminant function, from TYPE input
%       \'name1\'     name of first group of pair (group I)
%       \'name2\'     name of second group of pair (group J)
%       \'const\'     constant term of boundary equation (K)
%       \'linear\'    coefficients of linear term of boundary equation (L)
%       \'quadratic\' coefficient matrix of quadratic terms (Q)
%
%   For the \'linear\' and \'diaglinear\' types, the \'quadratic\' field is
%   absent, and a row x from the SAMPLE array is classified into group I
%   rather than group J if
%         0 < K + x*L
%   For the other types, x is classified into group I if
%         0 < K + x*L + x*Q*x\'
%
%   Example:
%      % Classify Fisher iris data using quadratic discriminant function
%      load fisheriris
%      x = meas(51:end,1:2);  % for illustrations use 2 species, 2 columns
%      y = species(51:end);
%      [c,err,post,logl,str] = classify(x,x,y,\'quadratic\');
%      h = gscatter(x(:,1),x(:,2),y,\'rb\',\'v^\');
%
%      % Classify a grid of values
%      [X,Y] = meshgrid(linspace(4.3,7.9), linspace(2,4.4));
%      X = X(:); Y = Y(:);
%      C = classify([X Y],x,y,\'quadratic\');
%      hold on; gscatter(X,Y,C,\'rb\',\'.\',1,\'off\'); hold off
%
%      % Draw boundary between two regions
%      hold on
%      K = str(1,2).const;
%      L = str(1,2).linear;
%      Q = str(1,2).quadratic;
%      % Plot the curve K + [x,y]*L + [x,y]*Q*[x,y]\' = 0:
%      f = @(x,y) K + L(1)*x + L(2)*y ...
%                   + Q(1,1)*x.^2 + (Q(1,2)+Q(2,1))*x.*y + Q(2,2)*y.^2;
%      ezplot(f,[4 8 2 4.5]);
%      hold off
%      title(\'Classification of Fisher iris data\')
%      legend(h)
%
%   See also fitcdiscr, fitcnb.

%   Copyright 1993-2018 The MathWorks, Inc.


%   References:
%     [1] Krzanowski, W.J., Principles of Multivariate Analysis,
%         Oxford University Press, Oxford, 1988.
%     [2] Seber, G.A.F., Multivariate Observations, Wiley, New York, 1984.

if nargin < 3
    error(message(\'stats:classify:TooFewInputs\'));
end

if nargin>3
    type = convertStringsToChars(type);
end

% grp2idx sorts a numeric grouping var ascending, and a string grouping
% var by order of first occurrence
[gindex,groups,glevels] = grp2idx(group);
nans = find(isnan(gindex));
if ~isempty(nans)
    training(nans,:) = [];
    gindex(nans) = [];
end
ngroups = length(groups);
gsize = hist(gindex,1:ngroups);
nonemptygroups = find(gsize>0);
nusedgroups = length(nonemptygroups);
if ngroups > nusedgroups
    warning(message(\'stats:classify:EmptyGroups\'));

end

[n,d] = size(training);
if size(gindex,1) ~= n
    error(message(\'stats:classify:TrGrpSizeMismatch\'));
elseif isempty(sample)
    sample = zeros(0,d,class(sample));  % accept any empty array but force correct size
elseif size(sample,2) ~= d
    error(message(\'stats:classify:SampleTrColSizeMismatch\'));
end
m = size(sample,1);

if nargin < 4 || isempty(type)
    type = \'linear\';
elseif ischar(type)
    types = {\'linear\',\'quadratic\',\'diaglinear\',\'diagquadratic\',\'mahalanobis\'};
    type = internal.stats.getParamVal(type,types,\'TYPE\');
else
    error(message(\'stats:classify:BadType\'));
end

% Default to a uniform prior
if nargin < 5 || isempty(prior)
    prior = ones(1, ngroups) / nusedgroups;
    prior(gsize==0) = 0;
    % Estimate prior from relative group sizes
elseif ischar(prior) && strncmpi(prior,\'empirical\',length(prior))
    %~isempty(strmatch(lower(prior), \'empirical\'))
    prior = gsize(:)\' / sum(gsize);
    % Explicit prior
elseif isnumeric(prior)
    if min(size(prior)) ~= 1 || max(size(prior)) ~= ngroups
        error(message(\'stats:classify:GrpPriorSizeMismatch\'));
    elseif any(prior < 0)
        error(message(\'stats:classify:BadPrior\'));
    end
    %drop empty groups
    prior(gsize==0)=0;
    prior = prior(:)\' / sum(prior); % force a normalized row vector
elseif isstruct(prior)
    [pgindex,pgroups] = grp2idx(prior.group);
   
    ord = NaN(1,ngroups);
    for i = 1:ngroups
      j = find(strcmp(groups(i), pgroups(pgindex)));
        if ~isempty(j)
            ord(i) = j;
        end
    end
    if any(isnan(ord))
        error(message(\'stats:classify:PriorBadGrpup\'));
    end
    prior = prior.prob(ord);
    if any(prior < 0)
        error(message(\'stats:classify:PriorBadProb\'));
    end
    prior(gsize==0)=0;
    prior = prior(:)\' / sum(prior); % force a normalized row vector
else
    error(message(\'stats:classify:BadPriorType\'));
end

% Add training data to sample for error rate estimation
if nargout > 1
    sample = [sample; training];
    mm = m+n;
else
    mm = m;
end

gmeans = NaN(ngroups, d);
for k = nonemptygroups
    gmeans(k,:) = mean(training(gindex==k,:),1);
end

D = NaN(mm, ngroups);
isquadratic = false;
switch type
    case \'linear\'
        if n <= nusedgroups
            error(message(\'stats:classify:NTrainingTooSmall\'));
        end
        % Pooled estimate of covariance.  Do not do pivoting, so that A can be
        % computed without unpermuting.  Instead use SVD to find rank of R.
        [Q,R] = qr(training - gmeans(gindex,:), 0);
        R = R / sqrt(n - nusedgroups); % SigmaHat = R\'*R
        s = svd(R);
        if any(s <= max(n,d) * eps(max(s)))
            error(message(\'stats:classify:BadLinearVar\'));
        end
        logDetSigma = 2*sum(log(s)); % avoid over/underflow
        
        % MVN relative log posterior density, by group, for each sample
        for k = nonemptygroups
            A = bsxfun(@minus,sample, gmeans(k,:)) / R;
            D(:,k) = log(prior(k)) - .5*(sum(A .* A, 2) + logDetSigma);
        end
        
    case \'diaglinear\'
        if n <= nusedgroups
            error(message(\'stats:classify:NTrainingTooSmall\'));
        end
        % Pooled estimate of variance: SigmaHat = diag(S.^2)
        S = std(training - gmeans(gindex,:)) * sqrt((n-1)./(n-nusedgroups));
        
        if any(S <= n * eps(max(S)))
            error(message(\'stats:classify:BadDiagLinearVar\'));
        end
        logDetSigma = 2*sum(log(S)); % avoid over/underflow
        
        if nargout >= 5
            R = S\';
        end
        
        % MVN relative log posterior density, by group, for each sample
        for k = nonemptygroups
            A=bsxfun(@times, bsxfun(@minus,sample,gmeans(k,:)),1./S);
            D(:,k) = log(prior(k)) - .5*(sum(A .* A, 2) + logDetSigma);
        end
        
    case {\'quadratic\' \'mahalanobis\'}
        if any(gsize == 1)
            error(message(\'stats:classify:BadTraining\'));
        end
        isquadratic = true;
        logDetSigma = zeros(ngroups,1,class(training));
        if nargout >= 5
          R = zeros(d,d,ngroups,class(training));
        end
        for k = nonemptygroups
            % Stratified estimate of covariance.  Do not do pivoting, so that A
            % can be computed without unpermuting.  Instead use SVD to find rank
            % of R.
            [Q,Rk] = qr(bsxfun(@minus,training(gindex==k,:),gmeans(k,:)), 0);
            Rk = Rk / sqrt(gsize(k) - 1); % SigmaHat = R\'*R
            s = svd(Rk);
%             if any(s <= max(gsize(k),d) * eps(max(s)))
%                 error(message(\'stats:classify:BadQuadVar\'));
%             end
            logDetSigma(k) = 2*sum(log(s)); % avoid over/underflow
            
            A = bsxfun(@minus, sample, gmeans(k,:)) /Rk;
            switch type
                case \'quadratic\'
                    % MVN relative log posterior density, by group, for each sample
                    D(:,k) = log(prior(k)) - .5*(sum(A .* A, 2) + logDetSigma(k));
                case \'mahalanobis\'
                    % Negative squared Mahalanobis distance, by group, for each
                    % sample.  Prior probabilities are not used
                    D(:,k) = -sum(A .* A, 2);
            end
            if nargout >=5 && ~isempty(Rk) 
                R(:,:,k) = inv(Rk);
            end
        end
        
    case \'diagquadratic\'
        if any(gsize == 1)
            error(message(\'stats:classify:BadTraining\'));
        end
        isquadratic = true;
        logDetSigma = zeros(ngroups,1,class(training));
        if nargout >= 5
            R = zeros(d,1,ngroups,class(training));
        end
        for k = nonemptygroups
            % Stratified estimate of variance:  SigmaHat = diag(S.^2)
            S = std(training(gindex==k,:));
            if any(S <= gsize(k) * eps(max(S)))
                error(message(\'stats:classify:BadDiagQuadVar\'));
            end
            logDetSigma(k) = 2*sum(log(S)); % avoid over/underflow
            
            % MVN relative log posterior density, by group, for each sample
            A=bsxfun(@times, bsxfun(@minus,sample,gmeans(k,:)),1./S);
            D(:,k) = log(prior(k)) - .5*(sum(A .* A, 2) + logDetSigma(k));
            if nargout >= 5
              R(:,:,k) = 1./S\';
            end
        end
end

% find nearest group to each observation in sample data
[maxD,outclass] = max(D, [], 2);

% Compute apparent error rate: percentage of training data that
% are misclassified, weighted by the prior probabilities for the groups.
if nargout > 1
    trclass = outclass(m+(1:n));
    outclass = outclass(1:m);
    
    miss = trclass ~= gindex;
    e = zeros(ngroups,1);
    for k = nonemptygroups
        e(k) = sum(miss(gindex==k)) / gsize(k);
    end
    err = prior*e;
end

if nargout > 2
    if strcmp(type, \'mahalanobis\')
        % Mahalanobis discrimination does not use the densities, so it\'s
        % possible that the posterior probs could disagree with the
        % classification.
        posterior = [];
        logp = [];
    else
        % Bayes\' rule: first compute p{x,G_j} = p{x|G_j}Pr{G_j} ...
        % (scaled by max(p{x,G_j}) to avoid over/underflow)
        P = exp(bsxfun(@minus,D(1:m,:),maxD(1:m)));
        sumP = nansum(P,2);
        % ... then Pr{G_j|x) = p(x,G_j} / sum(p(x,G_j}) ...
        % (numer and denom are both scaled, so it cancels out)
        posterior = bsxfun(@times,P,1./(sumP));
        if nargout > 3
            % ... and unconditional p(x) = sum(p(x,G_j}).
            % (remove the scale factor)
            logp = log(sumP) + maxD(1:m) - .5*d*log(2*pi);
        end
    end
end

%Convert outclass back to original grouping variable type
 outclass = glevels(outclass,:);

if nargout>=5
    pairs = combnk(nonemptygroups,2)\';
    npairs = size(pairs,2);
    K = zeros(1,npairs,class(training));
    L = zeros(d,npairs,class(training));
    if ~isquadratic
        % Methods with equal covariances across groups
        for j=1:npairs
            % Compute const (K) and linear (L) coefficients for
            % discriminating between each pair of groups
            i1 = pairs(1,j);
            i2 = pairs(2,j);
            mu1 = gmeans(i1,:)\';
            mu2 = gmeans(i2,:)\';
            if ~strcmp(type,\'diaglinear\')
                b = R \ ((R\') \ (mu1 - mu2));
            else
                b = (1./R).^2 .*(mu1 -mu2);
            end
            L(:,j) = b;
            K(j) = 0.5 * (mu1 + mu2)\' * b;
        end
    else
        % Methods with separate covariances for each group
        Q = zeros(d,d,npairs,class(training));
        for j=1:npairs
            % As above, but compute quadratic (Q) coefficients as well
            i1 = pairs(1,j);
            i2 = pairs(2,j);
            mu1 = gmeans(i1,:)\';
            mu2 = gmeans(i2,:)\';
            R1i = R(:,:,i1);    % note here the R array contains inverses
            R2i = R(:,:,i2);
            if ~strcmp(type,\'diagquadratic\')
                Rm1 = R1i\' * mu1;
                Rm2 = R2i\' * mu2;
            else
                Rm1 = R1i .* mu1;
                Rm2 = R2i .* mu2;
            end
            K(j) = 0.5 * (sum(Rm1.^2) - sum(Rm2.^2));
            if ~strcmp(type, \'mahalanobis\')
                K(j) = K(j) + 0.5 * (logDetSigma(i1)-logDetSigma(i2));
            end
            if ~strcmp(type,\'diagquadratic\')
                L(:,j) = (R1i*Rm1 - R2i*Rm2);
                Q(:,:,j) = -0.5 * (R1i*R1i\' - R2i*R2i\');
            else
                L(:,j) = (R1i.*Rm1 - R2i.*Rm2);
                Q(:,:,j) = -0.5 * diag(R1i.^2 - R2i.^2);
            end
        end
    end
    
    % For all except Mahalanobis, adjust for the priors
    if ~strcmp(type, \'mahalanobis\')
        K = K - log(prior(pairs(1,:))) + log(prior(pairs(2,:)));
    end
    
    % Return information as a structure
    coeffs = struct(\'type\',repmat({type},ngroups,ngroups));
    for k=1:npairs
        i = pairs(1,k);
        j = pairs(2,k);
        coeffs(i,j).name1 = glevels(i,:);
        coeffs(i,j).name2 = glevels(j,:);
        coeffs(i,j).const = -K(k);
        coeffs(i,j).linear = L(:,k);
        
        coeffs(j,i).name1 = glevels(j,:);
        coeffs(j,i).name2 = glevels(i,:);
        coeffs(j,i).const = K(k);
        coeffs(j,i).linear = -L(:,k);
        
        if isquadratic
            coeffs(i,j).quadratic = Q(:,:,k);
            coeffs(j,i).quadratic = -Q(:,:,k);
        end
    end
end

　　运行代码：

% 样本数据
training=[23,1,7,2,31,6.60,0.34,1.71
34,1,17,3,59,8.00,1.81,2.91
42,2,7,23,41,4.60,0.94,.94
39,1,19,5,48,13.10,1.93,4.36
35,1,9,1,34,5.00,0.40,1.30
37,1,1,3,24,15.10,1.80,1.82
29,1,13,1,42,7.40,1.46,1.65
32,2,11,6,75,23.30,7.76,9.72
28,2,2,3,23,6.40,0.19,1.29
26,1,4,3,27,10.50,2.47,.36];
% 用于构造判别函数的样本数据矩阵
group=[1,1,1,1,1,2,2,2,2,2]\';
% 参数group是与training相应的分组变量
sample=[53.00,1.00,9.00,18.00,50.00,11.20,2.02,3.58];
% 使用线性判别法分类
[class2,err2]=classify(sample,training,group,\'linear\')
% 使用马氏距离判别法分类
[class1,err1]=classify(sample,training,group,\'mahalanobis\')
% 使用二次判别法分类
[class3,err3]=classify(sample,training,group,\'quadratic\') 

　　运行示例：
>> exp64

class2 =

     1

err2 =

     0

警告: 秩亏，秩 = 4，tol =  1.864663e-14。

> In classify (line 286)

  In exp64 (line 19)

警告: 秩亏，秩 = 4，tol =  3.243458e-14。

> In classify (line 286)

  In exp64 (line 19)

class1 =

     1

err1 =

    0.2000

警告: 秩亏，秩 = 4，tol =  1.864663e-14。

> In classify (line 286)

  In exp64 (line 21)

警告: 秩亏，秩 = 4，tol =  3.243458e-14。

> In classify (line 286)

  In exp64 (line 21)

class3 =

     1

err3 =

    0.2000

　　显然，线性判别、二次判别和马氏距离判别都把样本归为第1类“已履行还贷责任”。那么最终结果该贷款人是已履行还贷责任一类。这里，因为在第二步MATLAB运行过程中，出现警告，这里本报告手写一份马氏距离判别的程序以作验证，二次判别以此类推。函数代码：

% 马氏距离判别
% 输入：G1组别1，G2组别2，待分类x
% 输出：分类序号
function w=ma(G1,G2,x)
[n1,m1]=size(G1);
[n2,m2]=size(G2);
if m1~=m2
    error(\'两个样本的列数应该相等！！\');
end
[n,~]=size(x);
w=zeros(n,1);
% 协方差
g1=cov(G1)./(n1-1);
g2=cov(G2)./(n2-1);
% 参数估计
for i=1:m1
    niu1(i)=mean(G1(:,i));
end
for i=1:m2
    niu2(i)=mean(G2(:,i));
end
% 马氏距离求解
v1=inv(g1);
v2=inv(g2);
d1=(x-niu1)\'.*v1.*(x-niu1);
d2=(x-niu2)\'.*v2.*(x-niu2);
k=1;
for j=1:n
    if det(d1)>det(d2)
        w(k)=2;
        k=k+1;
    else
        w(k)=1;
        k=k+1;
    end
end
end

　　运行代码：

G1=[23,1,7,2,31,6.60,0.34,1.71
34,1,17,3,59,8.00,1.81,2.91
42,2,7,23,41,4.60,0.94,.94
39,1,19,5,48,13.10,1.93,4.36
35,1,9,1,34,5.00,0.40,1.30];
G2=[37,1,1,3,24,15.10,1.80,1.82
29,1,13,1,42,7.40,1.46,1.65
32,2,11,6,75,23.30,7.76,9.72
28,2,2,3,23,6.40,0.19,1.29
26,1,4,3,27,10.50,2.47,.36];
x=[53.00,1.00,9.00,18.00,50.00,11.20,2.02,3.58];
w=ma(G1,G2,x)

　　运行示例：
>> exp64

警告: 矩阵接近奇异值，或者缩放错误。结果可能不准确。RCOND =  2.592596e-20。

> In ma (line 23)

  In exp64 (line 44)

警告: 矩阵接近奇异值，或者缩放错误。结果可能不准确。RCOND =  1.123714e-19。

> In ma (line 24)

  In exp64 (line 44)

w =

     1

　　明显，手写的马氏距离代码运行结果同样把样本归为第1类“已履行还贷责任”。但是，同样在手写的代码中同样有与classify相似的警告，从这次警告可以知道警告出现在协方差矩阵的逆接近奇异值。（并非是因为正定性的原因，之所以判断为非正定是因为值太小，被计算机认为等于0）

　　因此可以放大协方差矩阵的逆矩阵，已通过未修改的classify函数。这里不予详细展示。基于“一个半正定阵加一个正定阵就会成为一个正定阵”的思想，这里对手写的马氏距离判别代码作出如下修改：（原矩阵+一个1e-4的单位阵，当然可以更小，这里不予论证！）

% 马氏距离判别
% 输入：G1组别1，G2组别2，待分类x
% 输出：分类序号
function w=ma(G1,G2,x)
[n1,m1]=size(G1);
[n2,m2]=size(G2);
if m1~=m2
    error(\'两个样本的列数应该相等！！\');
end
[n,~]=size(x);
w=zeros(n,1);
% 协方差
g1=cov(G1)./(n1-1);
g2=cov(G2)./(n2-1);
[t1,t2]=size(g1);
g1=g1+eye(t1,t2)*1e-4;%修
g2=g2+eye(t1,t2)*1e-4;%修
% 参数估计
for i=1:m1
    niu1(i)=mean(G1(:,i));
end
for i=1:m2
    niu2(i)=mean(G2(:,i));
end
% 马氏距离求解
v1=inv(g1);
v2=inv(g2);
d1=(x-niu1)\'.*v1.*(x-niu1);
d2=(x-niu2)\'.*v2.*(x-niu2);
k=1;
for j=1:n
    if det(d1)>det(d2)
        w(k)=2;
        k=k+1;
    else
        w(k)=1;
        k=k+1;
    end
end
end

 

　　运行效果：

>> exp64

w =

1

　　如此，解决了奇异值的问题。那么是否可以把这种方法用在MATLAB自带的classify函数上呢？（笔者建议最好先做严格的数学论证，再行修改；笔者上述的修改是缺乏严格的数学论证的，待补！），这里不妨尝试一下，再修改一下MATLAB自带的classify函数。

 

 

 

 　　然后使用MATLAB求解问题2，代码如下

%例题2
clc,clear;
load(\'exp64.mat\');
x=exp64;
X=zscore(x);
std=corrcoef(X);
[vec,val]=eig(std);
newval=diag(val);
[y,i]=sort(newval);
rate=y/sum(y);
sumrate=0;
newi=[];
for k=length(y):-1:1
    sumrate=sumrate+rate(k);
    newi(length(y)+1-k)=i(k);
    if sumrate>0.85
        break;
    end
end
fprintf(\'主成分法：%g \n \n\',length(newi));
for i=1:1:length(newi)
    for j=1:1:length(y)
        aa(i,j)=sqrt(newval(newi(i)))*vec(j,newi(i));
    end
end
aaa=aa.*aa;
for i=1:1:length(newi)
    for j=1:1:length(y)
        zcfhz(i,j)=aa(i,j)/sqrt(sum(aaa(i,:)));
    end
end
fprintf(\'主成分载荷：\n\'),zcfhz

　　运行示例：

主成分法：1

主成分载荷：

zcfhz =

    0.4970    0.5146    0.4809    0.5069

解得这30名学生的第1主成分为：

 小结

　　在解第1题时，遇到了一些问题；得到的结论是并非是所有的样本数据都应该只用某种判别方法进行实验，有一些数据就不适合使用某种判别法，因此认识判别法的数学原理是必要的的。