import?numpy?as?np
import?pandas?as?pd
from?scipy?import?linalgstats
import?matplotlib.pyplot?as?plt
import?seaborn?as?sns
sns.set（）


#?產(chǎn)生訓(xùn)練數(shù)據(jù)
num_sample?=?100
a?=?10*np.random.randn（num_sample1）
x1?=?a+np.random.randn（num_sample1）
x2?=?1*np.sin（a）+np.random.randn（num_sample1）
x3?=?5*np.cos（5*a）+np.random.randn（num_sample1）
x4?=?0.8*x2+0.1*x3+np.random.randn（num_sample1）
x?=?np.hstack（（x1x2x3x4））
xx_train?=?x

#?產(chǎn)生測試數(shù)據(jù)
a?=?10*np.random.randn（num_sample1）
x1?=?a+np.random.randn（num_sample1）
x2?=?1*np.sin（a）+np.random.randn（num_sample1）
x3?=?5*np.cos（5*a）+np.random.randn（num_sample1）
x4?=?0.8*x2+0.1*x3+np.random.randn（num_sample1）
xx_test?=?np.hstack（（x1x2x3x4））
xx_test[50:1]?=?xx_test[50:1]+15*np.ones（50）

Xtrain?=?xx_train
Xtest?=?xx_test

#?標(biāo)準(zhǔn)化處理
X_mean?=?np.mean（Xtrain?axis=0）??#按列求Xtrain平均值
X_std?=?np.std（Xtrain?axis=0）??#求標(biāo)準(zhǔn)差
X_rowX_col?=?Xtrain.shape??#求Xtrain行、列數(shù)

Xtrain?=?（Xtrain-np.tile（X_mean（X_row1）））/np.tile（X_std（X_row1））

#?求協(xié)方差矩陣
sigmaXtrain?=?np.cov（Xtrain?rowvar=False）
#對協(xié)方差矩陣進(jìn)行特征分解，lamda為特征值構(gòu)成的對角陣，T的列為單位特征向量，且與lamda中的特征值一一對應(yīng)：
lamdaT?=?linalg.eig（sigmaXtrain）

#取對角元素（結(jié)果為一列向量），即lamda值，并上下反轉(zhuǎn)使其從大到小排列，主元個數(shù)初值為1，若累計貢獻(xiàn)率小于90%則增加主元個數(shù)
#?D?=?flipud（diag（lamda））
T?=?T[:lamda.argsort（）]??#將T的列按照lamda升序排列
lamda.sort（）??#將lamda按升序排列
D?=?-np.sort（-np.real（lamda））??#提取實數(shù)部分，按降序排列
num_pc?=?1
while?sum（D[0:num_pc]）/sum（D）?????num_pc?+=?1

#取與lamda相對應(yīng)的特征向量
P?=?T[:X_col-num_pc:X_col]
TT?=?Xtrain.dot（T）
TT1?=?Xtrain.dot（P）

#?求置信度為95%時的T2統(tǒng)計控制限
T2UCL1?=?num_pc*（X_row-1）*（X_row+1）*stats.f.ppf（0.95num_pcX_row-num_pc）/（X_row*（X_row-num_pc））

#?置信度為95%的Q統(tǒng)計控制限
theta?=?[]
for?i?in?range（1?num_pc+1）:
????theta.append（sum（（D[num_pc:]）**i））
h0?=?1?-?2*theta[0]*theta[2]/（3*theta[1]**2）
ca?=?stats.norm.ppf（0.9501）
QUCL?=?theta[0]*（h0*ca*np.sqrt（2*theta[1]）/theta[0]?+?1?+?theta[1]*h0*（h0?-?1）/theta[0]**2）**（1/h0）


#?在線監(jiān)測：
#?標(biāo)準(zhǔn)化處理
#?X_mean?=?np.mean（Xtest?