主成分分析方法中临界故障检测幅值计算代码

主成分分析方法中T^2统计量和SPE统计量的临界故障检测幅值计算的代码。

clear all
close all
clc

% * Read data
train_data = Creatmat（10000）;

% Suppose train_data has been read or created which has m rows and n
% columns.  Each row is a sample in which has n variables/sensors.
[m n] = size（train_data）;

% * Build PCA model with train_data
S = corrcoef（zscore（train_data））;

% Standardize train_data
pca_mean = mean（train_data）;
pca_std = std（train_data）;
[pc score latent] = princomp（zscore（train_data））;
sqrt_latent = sqrt（latent）;

f_t = zeros（1 n）;
f_q = zeros（1 n）;
for i = 1 : n 
    % P is the loading matrix P has the size of （n x i）
    P = pc（: 1 : i）;
    % tilde_P has the size of （n x （n - i））
    tilde_P = pc（: i + 1 : end）;

    % * Calculate the control limit of Q and T^2
    % delta_q and delta_t are control limits of Q and T^2 indices
    beta = 0.99;    
    delta_t = i * （m - 1） / （ m - i） * finv（beta i m-i）;
    theta = zeros（3 1）;
    for x = 1 : 3
        for y = i + 1 : n
            theta（x） = theta（x） + latent（y）^（x）;
        end
    end
    h = 1 - 2 * theta（1） * theta（3） / 3 / （theta（2）^2）;
    delta_q = theta（1） * （norminv（beta） * （2 * theta（2） * h^2）^0.5 / theta（1） + 1 + theta（2） * h * （h - 1） / theta（1）^2）^（1/h）;
    
    xi = eye（n）;
    M_M = diag（sqrt_latent（1 : i）） * P‘ * xi（: 1）;
    f_t（i） = 2 * （1 / （svd（M_M））） * delta_t;
    
    C_C = tilde_P * tilde_P‘ * xi（: 1）;
    f_q（i） = 2 * （1 / （svd（C_C））） * delta_q;
end
f_t
[min_f_t num_PC1] = min（f_t）
f_q
[min_f_q num_PC2] = min（f_q）

Y = Creatmat（1000）;
W = zscore（Y）;
[m2 n2] = size（Y）;

%******加入故障******                                  
b = 800;                                                                   %故障个数
% h = randn（b 1）;                                                         %加入干扰的类型：精度下降干扰
h = linspace（12 12 b）‘;                                                    %加入干扰的类型：幅值为1的恒值干扰
W（m2 - b + 1 : m2 1） = W（m2 - b + 1 : m2 1） + h;                         %叠加干扰
M = zscore（train_data）;
N = W;

%******建模******
S = cov（M）;                    
[coeff latent] = eig（S）;          
L = fliplr（diag（latent）‘）;     
U = fliplr（coeff）;
sqrtL = sqrt（L）;
V = sqrtL（ones（size（M 2） 1） :） * U;
%******PCA******
A1 = num_PC1;                                                                           
P1 = U（: 1 : A1）;                                                        %取前A列作为新特征向量

Ts1 = （diag（N * P1 * diag（1./ L（1 : A1）） * P1‘ * N‘））‘;                    %计算T2统计量
I1 = eye（size（P1 1））;
SPE1 = （diag（N * （I1 - P1 * P1‘） * （I1 - P1 * P1‘） * N‘））‘;                %计算SPE统计量

beta = 0.99;
TsBeta1 = A1 * （m - 1） / （m - A1） * finv（beta A1 m - A1）;                %求T2控制限                                                          
theta1 = zeros（3 1）;
for i = 1 : 3
    for j = A1 + 1 : size（M 2）
        theta1（i） = theta1（i） + L（j）^（i）;
    end
end
h1 = 1 - 2 * theta1（1） * theta1（3） / 3 / （theta1（2）^2）;
SPEBeta1 = theta1（1） * （norminv（

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     文件        7643  2014-07-17 21:22  A_LinJie_PCA.m
     文件         133  2014-01-29 18:02  Creatmat.m
     文件         252  2014-01-29 18:02  creat_mat.m