# How to implement the contribution analysis using PCA?

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
eigenvalues = eigenvalues

#### Working and correct:

# Contribution:
#### CDC:
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:
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:
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:
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):
# 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_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

PCA_scores = PCA_scores
eigenvalues = eigenvalues

# Equation 1:
Q_err = input_features - Q_eq_1
Q_res = np.sum(Q_err**2, axis=1)

# Equation 2:
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:
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:
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:
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:
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):
# # 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_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

PCA_scores = PCA_scores
eigenvalues = eigenvalues

Q_eq_2 = np.identity(7) - Q_eq_1
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):
# # 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_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)