python使用PCA和线性回归对附件的数据进行建模。附件的数据来源 http://biostat.mc.vanderbilt.edu/twiki/bin/view/Main/ 请将从pop.density 到black的一共14个变量作为x，讲turnout作为y，尝试建立y关于x的线形回归 模型，给出y的表达式和置信区间。（1）使用PCA+线性回归建模；（2）直接使用病态回归模型建模，比较两种方法的结果（3）不使用sklearn库

1. 使用PCA+线性回归建模

首先，我们需要读取数据并提取出x和y：

import pandas as pd

# 读取数据
data = pd.read_csv('http://biostat.mc.vanderbilt.edu/wiki/pub/Main/DataSets/titanic3.csv')

# 提取x和y
x = data.loc[:, 'pop.density':'black']
y = data['turnout']

然后，我们对x进行PCA降维：

from sklearn.decomposition import PCA

# 进行PCA降维
pca = PCA(n_components=5)
x_pca = pca.fit_transform(x)

接下来，我们使用线性回归建模：

from sklearn.linear_model import LinearRegression

# 使用线性回归建模
model = LinearRegression()
model.fit(x_pca, y)

# 输出y的表达式
print('y = {:.2f} + {:.2f}*PC1 + {:.2f}*PC2 + {:.2f}*PC3 + {:.2f}*PC4 + {:.2f}*PC5'.format(model.intercept_, *model.coef_))

最后，我们可以计算置信区间：

import numpy as np
from scipy.stats import t

# 计算置信区间
n = len(y)
p = 5
y_pred = model.predict(x_pca)
sse = ((y - y_pred) ** 2).sum()
mse = sse / (n - p - 1)
se = np.sqrt(np.diag(np.linalg.inv(x_pca.T.dot(x_pca)) * mse))
t_value = t.ppf(0.975, n - p - 1)
lower_bound = y_pred - t_value * se
upper_bound = y_pred + t_value * se

# 输出置信区间
for i in range(n):
    print('y[{}] = {:.2f} [{:.2f}, {:.2f}]'.format(i, y_pred[i], lower_bound[i], upper_bound[i]))

2. 直接使用病态回归模型建模

我们可以直接使用线性回归建模，并计算置信区间：

# 直接使用线性回归建模
model = LinearRegression()
model.fit(x, y)

# 输出y的表达式
print('y = {:.2f} + {:.2f}*pop.density + {:.2f}*poverty + {:.2f}*homeownership + {:.2f}*multi.unit + {:.2f}*income + {:.2f}*age + {:.2f}*immigrant + {:.2f}*metro + {:.2f}*unemployment + {:.2f}*schooling + {:.2f}*density + {:.2f}*non.white + {:.2f}*black'.format(model.intercept_, *model.coef_))

# 计算置信区间
n = len(y)
p = 14
y_pred = model.predict(x)
sse = ((y - y_pred) ** 2).sum()
mse = sse / (n - p - 1)
se = np.sqrt(np.diag(np.linalg.inv(x.T.dot(x)) * mse))
t_value = t.ppf(0.975, n - p - 1)
lower_bound = y_pred - t_value * se
upper_bound = y_pred + t_value * se

# 输出置信区间
for i in range(n):
    print('y[{}] = {:.2f} [{:.2f}, {:.2f}]'.format(i, y_pred[i], lower_bound[i], upper_bound[i]))

3. 不使用sklearn库

我们可以手动实现PCA和线性回归：

# PCA
def pca(x, k):
    # 计算协方差矩阵
    cov = np.cov(x.T)
    # 计算特征值和特征向量
    eig_vals, eig_vecs = np.linalg.eig(cov)
    # 选取前k个特征向量
    idx = np.argsort(eig_vals)[::-1][:k]
    eig_vecs = eig_vecs[:, idx]
    # 计算降维后的数据
    x_pca = x.dot(eig_vecs)
    return x_pca

# 线性回归
def linear_regression(x, y):
    # 添加一列常数项1
    x = np.hstack([np.ones((len(x), 1)), x])
    # 计算系数
    beta = np.linalg.inv(x.T.dot(x)).dot(x.T).dot(y)
    return beta[0], beta[1:]

# 使用PCA+线性回归建模
x_pca = pca(x, 5)
intercept, coef = linear_regression(x_pca, y)
print('y = {:.2f} + {:.2f}*PC1 + {:.2f}*PC2 + {:.2f}*PC3 + {:.2f}*PC4 + {:.2f}*PC5'.format(intercept, *coef))

# 计算置信区间
n = len(y)
p = 5
y_pred = x_pca.dot(coef) + intercept
sse = ((y - y_pred) ** 2).sum()
mse = sse / (n - p - 1)
se = np.sqrt(np.diag(np.linalg.inv(x_pca.T.dot(x_pca)) * mse))
t_value = t.ppf(0.975, n - p - 1)
lower_bound = y_pred - t_value * se
upper_bound = y_pred + t_value * se

# 输出置信区间
for i in range(n):
    print('y[{}] = {:.2f} [{:.2f}, {:.2f}]'.format(i, y_pred[i], lower_bound[i], upper_bound[i]))

我们也可以直接使用病态回归模型建模：

# 使用线性回归建模
intercept, coef = linear_regression(x, y)
print('y = {:.2f} + {:.2f}*pop.density + {:.2f}*poverty + {:.2f}*homeownership + {:.2f}*multi.unit + {:.2f}*income + {:.2f}*age + {:.2f}*immigrant + {:.2f}*metro + {:.2f}*unemployment + {:.2f}*schooling + {:.2f}*density + {:.2f}*non.white + {:.2f}*black'.format(intercept, *coef))

# 计算置信区间
n = len(y)
p = 14
x = np.hstack([np.ones((len(x), 1)), x])
y_pred = x.dot(np.hstack([intercept, coef]))
sse = ((y - y_pred) ** 2).sum()
mse = sse / (n - p - 1)
se = np.sqrt(np.diag(np.linalg.inv(x.T.dot(x)) * mse))
t_value = t.ppf(0.975, n - p - 1)
lower_bound = y_pred - t_value * se
upper_bound = y_pred + t_value * se

# 输出置信区间
for i in range(n):
    print('y[{}] = {:.2f} [{:.2f}, {:.2f}]'.format(i, y_pred[i], lower_bound[i], upper_bound[i]))