0
$\begingroup$

I have been looking into implementing the Q-Residual and Hotelling's T statistics calculation to the PCA components which is similar to the following article and website:

So, I tried to create the following code where it implements the calculation, however, I am not feeling confident around converting the equation into code, wherein CDC calculation part I keep on getting zeros across all Q, T and combined. What is my mistake in my code?

Side Question:

Will it be possible to implement the contribution analysis using LDA (Linear Discriminant Analysis)?

Code:

def hotelling_tsquared_index_PCA(input_features, PCA_scores, PCA_loading, eigenvalues):
    PCA_scores = PCA_scores
    PCA_loading = PCA_loading
    eigenvalues = eigenvalues

    #### Working and correct:    
    # t2 = np.array(np.diag(np.linalg.multi_dot([input_features, PCA_loading, np.linalg.inv(np.diag(eigenvalues)), PCA_loading.T, input_features.T])))
    t2 = np.array([xi.dot(PCA_loading).dot(np.diag(eigenvalues**-1)).dot(PCA_loading.T).dot(xi.T) for xi in input_features])

    # Contribution:
    #### CDC:
    # t2_CDC_eq_1 = np.linalg.multi_dot(PCA_loading, np.linalg.inv(np.diag(eigenvalues)), PCA_loading.T)
    t2_CDC_eq_1 = np.linalg.multi_dot([PCA_loading, np.linalg.inv(np.diag(eigenvalues)), PCA_loading.T])
    eye_M = np.identity(7)
    t2_CDC = np.array(np.diag(np.linalg.multi_dot([input_features, (t2_CDC_eq_1 ** 0.5), eye_M.T, eye_M, (t2_CDC_eq_1 ** 0.5), input_features.T])))
    t2_CDC = np.where(np.isnan(t2_CDC), 0, t2_CDC)
    
    #### PDC:
    # t2_PDC_eq_1 = np.linalg.multi_dot(PCA_loading, np.linalg.inv(np.diag(eigenvalues)), PCA_loading.T)
    t2_PDC_eq_1 = np.linalg.multi_dot([PCA_loading, np.linalg.inv(np.diag(eigenvalues)), PCA_loading.T])
    eye_M = np.identity(7)
    t2_PDC = np.array(np.diag(np.linalg.multi_dot([input_features, t2_PDC_eq_1, eye_M.T, eye_M, input_features.T])))
    t2_PDC = np.where(np.isnan(t2_PDC), 0, t2_PDC)
    
    #### DC:
    # t2_DC_eq_1 = np.linalg.multi_dot(PCA_loading, np.linalg.inv(np.diag(eigenvalues)), PCA_loading.T)
    t2_DC_eq_1 = np.linalg.multi_dot([PCA_loading, np.linalg.inv(np.diag(eigenvalues)), PCA_loading.T])
    eye_M = np.identity(7)
    t2_DC = np.array(np.diag(np.linalg.multi_dot([input_features, eye_M.T, eye_M, t2_DC_eq_1, eye_M.T, eye_M, input_features.T])))
    t2_DC = np.where(np.isnan(t2_DC), 0, t2_DC)
    
    #### RBC:
    # t2_RBC_eq_1 = np.linalg.multi_dot(PCA_loading, np.linalg.inv(np.diag(eigenvalues)), PCA_loading.T)
    t2_RBC_eq_1 = np.linalg.multi_dot([PCA_loading, np.linalg.inv(np.diag(eigenvalues)), PCA_loading.T])
    eye_M = np.identity(7)
    t2_RBC = np.array(np.diag(np.linalg.multi_dot([input_features, t2_RBC_eq_1, eye_M.T, np.linalg.pinv(np.linalg.multi_dot([eye_M, t2_RBC_eq_1, eye_M.T])), eye_M, t2_RBC_eq_1, input_features.T])))
    t2_RBC = np.where(np.isnan(t2_RBC), 0, t2_RBC)
    
    #### ABC (1):
    # # t2_ABC_eq_1 = np.linalg.multi_dot(PCA_loading, np.linalg.inv(np.diag(eigenvalues)), PCA_loading.T)
    # t2_ABC_eq_1 = np.linalg.multi_dot([PCA_loading, np.linalg.inv(np.diag(eigenvalues)), PCA_loading.T])
    # print(t2_ABC_eq_1)
    # print('Contribution Equation [1] of PCA hotellings T^2 matrix shape:', np.array(t2_ABC_eq_1).shape)
    # t2_ABC = t2_RBC / t2_ABC_eq_1
    # # print(t2_ABC)
    # print('Contribution of PCA hotellings T^2 matrix shape:', np.array(t2_ABC).shape)
    # t2_ABC = np.where(np.isnan(t2_ABC), 0, t2_ABC)
    # print(t2_ABC)

    # #### ABC (2):    
    # eye_M = np.identity(7)
    # t2_ABC_eq_1 = np.linalg.multi_dot([PCA_loading, np.linalg.inv(np.diag(eigenvalues)), input_features])
    # t2_ABC_eq_2 = (np.linalg.multi_dot([eye_M.T, t2_ABC_eq_1, PCA_loading.T]) ** 2)
    # t2_ABC_eq_3 = np.linalg.multi_dot([eye_M.T, t2_ABC_eq_1, eye_M, input_features.T, t2_ABC_eq_1, input_features])
    # print(t2_ABC_eq_1)
    # print('Contribution Equation [1] of PCA hotellings T^2 matrix shape:', np.array(t2_ABC_eq_1).shape)
    # t2_ABC = t2_ABC_eq_2 / t2_ABC_eq_3
    # # print(t2_ABC)
    # print('Contribution of PCA hotellings T^2 matrix shape:', np.array(t2_ABC).shape)
    # t2_ABC = np.where(np.isnan(t2_ABC), 0, t2_ABC)
    # print(t2_ABC)
    # return t2
    return t2, t2_CDC, t2_PDC, t2_DC, t2_RBC

def Q_Residual_index_PCA(input_features, PCA_scores, PCA_loading, eigenvalues):
    PCA_scores = PCA_scores
    PCA_loading = PCA_loading
    eigenvalues = eigenvalues
    
    # Equation 1:
    Q_eq_1 = np.dot(PCA_scores, PCA_loading.T)
    Q_err = input_features - Q_eq_1
    Q_res = np.sum(Q_err**2, axis=1)

    # Equation 2:
    Q_eq_1 = np.linalg.multi_dot([LDA_loading.T, LDA_loading])
    Q_eq_2 = np.identity(7) - Q_eq_1
    Q_res = np.linalg.multi_dot([input_features, Q_eq_2, input_features.T])
    
    # Contribution:
    #### CDC:
    # Q_res_CDC_eq_1 = np.linalg.multi_dot(PCA_loading, np.linalg.inv(np.diag(eigenvalues)), PCA_loading.T)
    Q_res_CDC_eq_1 = np.linalg.multi_dot([PCA_loading, np.linalg.inv(np.diag(eigenvalues)), PCA_loading.T])
    eye_M = np.identity(7)
    Q_res_CDC = np.array(np.diag(np.linalg.multi_dot([input_features, (Q_res_CDC_eq_1 ** 0.5), eye_M.T, eye_M, (Q_res_CDC_eq_1 ** 0.5), input_features.T])))
    Q_res_CDC = np.where(np.isnan(Q_res_CDC), 0, Q_res_CDC)
    
    #### PDC:
    # Q_res_PDC_eq_1 = np.linalg.multi_dot(PCA_loading, np.linalg.inv(np.diag(eigenvalues)), PCA_loading.T)
    Q_res_PDC_eq_1 = np.linalg.multi_dot([PCA_loading, np.linalg.inv(np.diag(eigenvalues)), PCA_loading.T])
    eye_M = np.identity(7)
    Q_res_PDC = np.array(np.diag(np.linalg.multi_dot([input_features, Q_res_PDC_eq_1, eye_M.T, eye_M, input_features.T])))
    Q_res_PDC = np.where(np.isnan(Q_res_PDC), 0, Q_res_PDC)
    
    #### DC:
    # Q_res_DC_eq_1 = np.linalg.multi_dot(PCA_loading, np.linalg.inv(np.diag(eigenvalues)), PCA_loading.T)
    Q_res_DC_eq_1 = np.linalg.multi_dot([PCA_loading, np.linalg.inv(np.diag(eigenvalues)), PCA_loading.T])
    eye_M = np.identity(7)
    Q_res_DC = np.array(np.diag(np.linalg.multi_dot([input_features, eye_M.T, eye_M, Q_res_DC_eq_1, eye_M.T, eye_M, input_features.T])))
    Q_res_DC = np.where(np.isnan(Q_res_DC), 0, Q_res_DC)
    
    #### RBC:
    # Q_res_RBC_eq_1 = np.linalg.multi_dot(PCA_loading, np.linalg.inv(np.diag(eigenvalues)), PCA_loading.T)
    Q_res_RBC_eq_1 = np.linalg.multi_dot([PCA_loading, np.linalg.inv(np.diag(eigenvalues)), PCA_loading.T])
    eye_M = np.identity(7)
    Q_res_RBC_eq_2 = np.linalg.multi_dot([eye_M, Q_res_RBC_eq_1, eye_M.T])
    Q_res_RBC = np.array(np.diag(np.linalg.multi_dot([input_features, Q_res_RBC_eq_1, eye_M.T, np.linalg.pinv(Q_res_RBC_eq_2), eye_M, Q_res_RBC_eq_1, input_features.T])))
    Q_res_RBC = np.where(np.isnan(Q_res_RBC), 0, Q_res_RBC)
    
    # #### ABC (1):
    # # # Q_res_ABC_eq_1 = np.linalg.multi_dot(PCA_loading, np.linalg.inv(np.diag(eigenvalues)), PCA_loading.T)
    # # Q_res_ABC_eq_1 = np.linalg.multi_dot([PCA_loading, np.linalg.inv(np.diag(eigenvalues)), PCA_loading.T])
    # # print(Q_res_ABC_eq_1)
    # # print('Contribution Equation [1] of PCA hotellings T^2 matrix shape:', np.array(Q_res_ABC_eq_1).shape)
    # # Q_res_ABC = Q_res_RBC / Q_res_ABC_eq_1
    # # # print(Q_res_ABC)
    # # print('Contribution of PCA hotellings T^2 matrix shape:', np.array(Q_res_ABC).shape)
    # # Q_res_ABC = np.where(np.isnan(Q_res_ABC), 0, Q_res_ABC)
    # # print(Q_res_ABC)

    # #### ABC (2):    
    # eye_M = np.identity(7)
    # Q_res_ABC_eq_1 = np.linalg.multi_dot([PCA_loading, np.linalg.inv(np.diag(eigenvalues)), input_features])
    # Q_res_ABC_eq_2 = (np.linalg.multi_dot([eye_M.T, Q_res_ABC_eq_1, PCA_loading.T]) ** 2)
    # Q_res_ABC_eq_3 = np.linalg.multi_dot([eye_M.T, Q_res_ABC_eq_1, eye_M, input_features.T, Q_res_ABC_eq_1, input_features])
    # print(Q_res_ABC_eq_1)
    # print('Contribution Equation [1] of PCA hotellings T^2 matrix shape:', np.array(Q_res_ABC_eq_1).shape)
    # Q_res_ABC = Q_res_ABC_eq_2 / Q_res_ABC_eq_3
    # # print(Q_res_ABC)
    # print('Contribution of PCA hotellings T^2 matrix shape:', np.array(Q_res_ABC).shape)
    # Q_res_ABC = np.where(np.isnan(Q_res_ABC), 0, Q_res_ABC)
    # print(Q_res_ABC)    
    return Q_res, Q_res_CDC, Q_res_PDC, Q_res_DC, Q_res_RBC

def Phi_index_PCA(input_features, PCA_scores, PCA_loading, eigenvalues):
    PCA_scores = PCA_scores
    PCA_loading = PCA_loading
    eigenvalues = eigenvalues
    
    Q_eq_1 = np.linalg.multi_dot([PCA_loading, PCA_loading.T])
    Q_eq_2 = np.identity(7) - Q_eq_1
    t2_eq_1 = np.linalg.multi_dot([PCA_loading, np.linalg.inv(np.diag(eigenvalues)), PCA_loading.T])
    Phi_eq_1 = Q_eq_2 + t2_eq_1
    Phi_res = np.array(np.diag(np.linalg.multi_dot([input_features, Phi_eq_1, input_features.T])))
    
    # Contribution:
    #### CDC:
    Phi_res_CDC_eq_1 = Phi_eq_1
    eye_M = np.identity(7)
    Phi_res_CDC = np.array(np.diag(np.linalg.multi_dot([input_features, (Phi_res_CDC_eq_1 ** 0.5), eye_M.T, eye_M, (Phi_res_CDC_eq_1 ** 0.5), input_features.T])))
    Phi_res_CDC = np.where(np.isnan(Phi_res_CDC), 0, Phi_res_CDC)
    
    #### PDC:
    Phi_res_PDC_eq_1 = Phi_eq_1
    eye_M = np.identity(7)
    Phi_res_PDC = np.array(np.diag(np.linalg.multi_dot([input_features, Phi_res_PDC_eq_1, eye_M.T, eye_M, input_features.T])))
    Phi_res_PDC = np.where(np.isnan(Phi_res_PDC), 0, Phi_res_PDC)
    
    #### DC:
    Phi_res_DC_eq_1 = Phi_eq_1
    eye_M = np.identity(7)
    Phi_res_DC = np.array(np.diag(np.linalg.multi_dot([input_features, eye_M.T, eye_M, Phi_res_DC_eq_1, eye_M.T, eye_M, input_features.T])))
    Phi_res_DC = np.where(np.isnan(Phi_res_DC), 0, Phi_res_DC)
    
    #### RBC:
    Phi_res_RBC_eq_1 = Phi_eq_1
    eye_M = np.identity(7)
    Phi_res_RBC = np.array(np.diag(np.linalg.multi_dot([input_features, Phi_res_RBC_eq_1, eye_M.T, np.linalg.pinv(np.linalg.multi_dot([eye_M, Phi_res_RBC_eq_1, eye_M.T])), eye_M, Phi_res_RBC_eq_1, input_features.T])))
    Phi_res_RBC = np.where(np.isnan(Phi_res_RBC), 0, Phi_res_RBC)
    
    # #### ABC (1):
    # # # Phi_res_ABC_eq_1 = np.linalg.multi_dot(PCA_loading, np.linalg.inv(np.diag(eigenvalues)), PCA_loading.T)
    # # Phi_res_ABC_eq_1 = np.linalg.multi_dot([PCA_loading, np.linalg.inv(np.diag(eigenvalues)), PCA_loading.T])
    # # print(Phi_res_ABC_eq_1)
    # # print('Contribution Equation [1] of PCA hotellings T^2 matrix shape:', np.array(Phi_res_ABC_eq_1).shape)
    # # Phi_res_ABC = Phi_res_RBC / Phi_res_ABC_eq_1
    # # # print(Phi_res_ABC)
    # # print('Contribution of PCA hotellings T^2 matrix shape:', np.array(Phi_res_ABC).shape)
    # # Phi_res_ABC = np.where(np.isnan(Phi_res_ABC), 0, Phi_res_ABC)
    # # print(Phi_res_ABC)

    # #### ABC (2):    
    # eye_M = np.identity(7)
    # Phi_res_ABC_eq_1 = np.linalg.multi_dot([PCA_loading, np.linalg.inv(np.diag(eigenvalues)), input_features])
    # Phi_res_ABC_eq_2 = (np.linalg.multi_dot([eye_M.T, Phi_res_ABC_eq_1, PCA_loading.T]) ** 2)
    # Phi_res_ABC_eq_3 = np.linalg.multi_dot([eye_M.T, Phi_res_ABC_eq_1, eye_M, input_features.T, Phi_res_ABC_eq_1, input_features])
    # print(Phi_res_ABC_eq_1)
    # print('Contribution Equation [1] of PCA hotellings T^2 matrix shape:', np.array(Phi_res_ABC_eq_1).shape)
    # Phi_res_ABC = Phi_res_ABC_eq_2 / Phi_res_ABC_eq_3
    # # print(Phi_res_ABC)
    # print('Contribution of PCA hotellings T^2 matrix shape:', np.array(Phi_res_ABC).shape)
    # Phi_res_ABC = np.where(np.isnan(Phi_res_ABC), 0, Phi_res_ABC)
    # print(Phi_res_ABC)    
    return Phi_res, Phi_res_CDC, Phi_res_PDC, Phi_res_DC, Phi_res_RBC 

#### Part of Main Program: 
print('Model: ', pca)
# pca = PCA(n_components=2)
PrincipalComponents = pca.fit_transform(X_std)
pca_f = pca.fit(X_std)
PCA_scores = pca_f.transform(X_std)
PCA_loading = (pca_f.components_).T
eigenvalues = pca_f.explained_variance_
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.