# 必要な関連ライブラリをインポートする
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.cross_decomposition import PLSRegression
from sklearn.metrics import r2_score, mean_absolute_error
from sklearn.model_selection import cross_val_score

# データ準備および正規化の変換処理
# ----------------------------------------------------
# チェックしたい構造化データ読み込みと説明変数Xと目的変数yの生成
df_ref = pd.read_csv('test_data.csv')
y = df_ref['目標値'].to_numpy()
X = df_ref.loc[:,'説明変数A':'説明変数E'].to_numpy()

# 説明変数Xに対して、z-score正規化によるデータ標準化
X_mean = np.mean(X, axis=0)
X_std = np.std(X, axis=0)
X_zscore_normalized = (X - X_mean) / X_std

# 目標変数に対して、オリジナル値をそのまま利用する
y_zscore_normalized = y


# クロスバリデーションの検証
#--------------------------------------------------
# コンポーネント数の候補
n_components_list = [1, 2, 3, 4, 5]

# 平均二乗誤差（MSE）を格納するリスト
mse_list = []

for n in n_components_list:
    pls = PLSRegression(n_components=n)
    mse = -np.mean(cross_val_score(pls, X_zscore_normalized, y_zscore_normalized, cv=5, scoring='neg_mean_squared_error'))
    mse_list.append(mse)
    print(f'n_components = {n}, MSE = {mse:.4f}')

# 最適なコンポーネント数の選択
optimal_n = n_components_list[np.argmin(mse_list)]
print(f'最適なコンポーネント数: {optimal_n}')


# PLS回帰モデルの構築と学習
#--------------------------------------------------
pls = PLSRegression(n_components=optimal_n)
pls.fit(X_zscore_normalized, y_zscore_normalized)
y_pred = pls.predict(X_zscore_normalized).ravel()

r2 = r2_score(y_zscore_normalized, y_pred)
mae = mean_absolute_error(y_zscore_normalized, y_pred)
print(f"決定係数 R^2: {r2:.4f}")
print(f"平均絶対誤差 MAE: {mae:.4f}")

# 残差（予測値と実際値の偏差）の計算
residuals = y_zscore_normalized - y_pred

# 残差の平均二乗誤差（MSE）の計算
MSE = np.mean(residuals ** 2)

# スコア行列の取得
T = pls.x_scores_

# レバレッジの計算
leverages = np.sum(T ** 2, axis=1) / np.sum(T ** 2)

# 外部スチューデント化残差の計算
#--------------------------------------------------
n = n_samples
p = pls.n_components # 使用したコンポーネント数
RSS = np.sum(residuals ** 2)

externally_studentized_residuals = np.zeros(n)

for i in range(n):
    # i番目の観測値を除いた残差平方和と自由度
    RSS_i = RSS - residuals[i] ** 2
    #自由度参考資料
    # https://ja.wikipedia.org/wiki/%E3%82%B9%E3%83%81%E3%83%A5%E3%83%BC%E3%83%87%E3%83%B3%E3%83%88%E5%8C%96%E6%AE%8B%E5%B7%AE
    # https://www.business-research-lab.com/211105-2/
    df_i = n - p - 2 # i番目を除いたため、自由度が1減る
    s_i = np.sqrt(RSS_i / df_i)
    externally_studentized_residuals[i] = residuals[i] / (s_i * np.sqrt(1 - leverages[i]))

# プロットの作成
plt.figure(figsize=(10, 6))
plt.scatter(leverages, externally_studentized_residuals)
plt.xlabel('Leverage')
plt.ylabel('Externally Studentized Residual')
plt.title('Externally Studentized Residuals vs Leverage')
plt.grid(True)
plt.show()