import pandas as pd import numpy as np from sklearn.decomposition import PCA from sklearn.preprocessing import StandardScaler import matplotlib.pyplot as plt # 读取数据 data = pd.read_csv('D:/pythonProject/venv/BostonHousing2.csv') # 提取前13个指标的数据 X = data.iloc[:, 5:18].values # 数据标准化 scaler = StandardScaler() X_scaled = scaler.fit_transform(X) # 主成分分析 pca = PCA() X_pca = pca.fit_transform(X_scaled) # 特征值和特征向量 eigenvalues = pca.explained_variance_ eigenvectors = pca.components_.T # 碎石图 variance_explained = np.cumsum(eigenvalues / np.sum(eigenvalues)) plt.plot(range(6, 19), variance_explained, marker='o') plt.xlabel('Number of Components') plt.ylabel('Cumulative Proportion of Variance Explained') plt.title('Scree Plot') plt.show() # 选择主成分个数 n_components = np.sum(variance_explained <= 0.95) + 1 # 前2个主成分的载荷图 loadings = pd.DataFrame(eigenvectors[:, 0:2], columns=['PC1', 'PC2'], index=data.columns[0:13]) plt.figure(figsize=(10, 6)) plt.scatter(loadings['PC1'], loadings['PC2'], alpha=0.7) for i, feature in enumerate(loadings.index): plt.text(loadings['PC1'][i], loadings['PC2'][i], feature) plt.xlabel('PC1') plt.ylabel('PC2') plt.title('Loading Plot') plt.grid() plt.show() # 主成分得分图 scores = pd.DataFrame(X_pca[:, 0:n_components], columns=['PC{}'.format(i+1) for i in range(n_components)]) plt.figure(figsize=(10, 6)) plt.scatter(scores['PC1'], scores['PC2'], alpha=0.7) for i, label in enumerate(data['MEDV']): plt.text(scores['PC1'][i], scores['PC2'][i], label) plt.xlabel('PC1') plt.ylabel('PC2') plt.title('Scores Plot') plt.grid() plt.show() # 综合评估和排序 data['PC1_score'] = X_pca[:, 0] sorted_data = data.sort_values(by='PC1_score') # 主成分回归模型 from sklearn.linear_model import LinearRegression Y = data['MEDV'].values.reshape(-1, 1) X_pca_regression = X_pca[:, 0].reshape(-1, 1) regression_model = LinearRegression() regression_model.fit(X_pca_regression, Y) # 回归方程 intercept = regression_model.intercept_[0] slope = regression_model.coef_[0][0] equation = "MEDV = {:.2f} + {:.2f} * PC1".format(intercept, slope) print("Regression Equation:", equation) # 最小二乘估计结果 from statsmodels.api import OLS X_const = np.concatenate((np.ones((506, 1)), X_pca_regression), axis=1) ols_model = OLS(Y, X_const).fit() print("OLS Regression Summary:") print(ols_model.summary())