判別分析與主成分分析

實驗目的

　　（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題時，遇到了一些問題；得到的結論是並非是所有的樣本數據都應該只用某種判別方法進行實驗，有一些數據就不適合使用某種判別法，因此認識判別法的數學原理是必要的的。