单变量线性回归

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

path =  'linear_regression.txt'
data = pd.read_csv(path, header=None, names=['Population', 'Profit'])
data.head()

	Population	Profit
0	6.1101	17.5920
1	5.5277	9.1302
2	8.5186	13.6620
3	7.0032	11.8540
4	5.8598	6.8233

data['Population'].describe()

count    97.000000
mean      8.159800
std       3.869884
min       5.026900
25%       5.707700
50%       6.589400
75%       8.578100
max      22.203000
Name: Population, dtype: float64

看下数据长什么样子

data.plot(kind='scatter', x='Population', y='Profit', figsize=(12,8))
#plt.figure(figsize=(12,8))
#plt.scatter(data['Population'],data['Profit'])
#lt.xlabel('Population')
#plt.ylabel('Porfit')
plt.show()

现在让我们使用梯度下降来实现线性回归，以最小化代价函数。

首先，我们将创建一个以参数θ为特征函数的代价函数 $$J\left( \theta \right)=\frac{1}{2m}\sum\limits_{i=1}^{m}{{{\left( {{h}_{\theta }}\left( {{x}^{(i)}} \right)-{{y}^{(i)}} \right)}^{2}}}$$ 其中 $${{h}_{\theta }}\left( x \right)={{\theta }^{T}}X={{\theta }_{0}}{{x}_{0}}+{{\theta }_{1}}{{x}_{1}}+{{\theta }_{2}}{{x}_{2}}+...+{{\theta }_{n}}{{x}_{n}}$$

def computeCost(X, y, theta):
    inner = np.power(((X * theta.T) - y), 2)
    return np.sum(inner) / (2 * len(X))

让我们在训练集中添加一列，以便可以使用向量化的解决方案来计算代价和梯度。

data.insert(0, 'Ones', 1)

data.head()

	Ones	Population	Profit
0	1	6.1101	17.5920
1	1	5.5277	9.1302
2	1	8.5186	13.6620
3	1	7.0032	11.8540
4	1	5.8598	6.8233

现在我们来做一些变量初始化。

# set X (training data) and y (target variable)
cols = data.shape[1]
X = data.iloc[:,0:cols-1]#X是所有行，去掉最后一列
y = data.iloc[:,cols-1:cols]#X是所有行，最后一列

观察下 X (训练集) and y (目标变量)是否正确.

X.head()#head()是观察前5行

	Ones	Population
0	1	6.1101
1	1	5.5277
2	1	8.5186
3	1	7.0032
4	1	5.8598

y.head()

	Profit
0	17.5920
1	9.1302
2	13.6620
3	11.8540
4	6.8233

代价函数是应该是numpy矩阵，所以我们需要转换X和Y，然后才能使用它们。 我们还需要初始化theta。

X = np.matrix(X.values)
y = np.matrix(y.values)
theta = np.matrix(np.array([0,0]))

theta 是一个(1,2)矩阵

theta

matrix([[0, 0]])

看下维度

X.shape, theta.shape, y.shape

((97, 2), (1, 2), (97, 1))

计算代价函数 (theta初始值为0).

computeCost(X, y, theta)

32.072733877455676

令截距为0时，画出代价函数随theta_1的变化

theta_1_list=np.linspace(-0.5,2.5,50)
cost_list=[]
for i in range(len(theta_1_list)):
    cost_list.append(computeCost(X,y,np.matrix(np.array([0,theta_1_list[i]]))))
plt.figure(figsize=(12,8))
plt.plot(theta_1_list,cost_list)
plt.xlabel('theta_1')
plt.ylabel('Cost(theta_1)')
plt.axvline(x=0,ls="--",c="red")
plt.show()

batch gradient decent（批量梯度下降） $${{\theta }_{j}}:={{\theta }_{j}}-\alpha \frac{\partial }{\partial {{\theta }_{j}}}J\left( \theta \right)$$

def gradientDescent(X, y, theta, alpha, iters):
    temp = np.matrix(np.zeros(theta.shape))
    parameters = int(theta.ravel().shape[1])
    cost = np.zeros(iters)
    
    for i in range(iters):
        error = (X * theta.T) - y
        
        for j in range(parameters):
            term = np.multiply(error, X[:,j])
            temp[0,j] = theta[0,j] - ((alpha / len(X)) * np.sum(term))
            
        theta = temp
        cost[i] = computeCost(X, y, theta)
        
    return theta, cost

初始化一些附加变量 - 学习速率α和要执行的迭代次数。

alpha = 0.01
iters = 1500

现在让我们运行梯度下降算法来将我们的参数θ适合于训练集。

g, cost = gradientDescent(X, y, theta, alpha, iters)
g

matrix([[-3.63029144,  1.16636235]])

最后，我们可以使用我们拟合的参数计算训练模型的代价函数（误差）。

computeCost(X, y, g)

4.483388256587726

现在我们来绘制线性模型以及数据，直观地看出它的拟合。

x = np.linspace(data.Population.min(), data.Population.max(), 100)
f = g[0, 0] + (g[0, 1] * x)

fig, ax = plt.subplots(figsize=(12,8))
ax.plot(x, f, 'r', label='Prediction')
ax.scatter(data.Population, data.Profit, label='Traning Data')
ax.legend(loc=2)
ax.set_xlabel('Population')
ax.set_ylabel('Profit')
ax.set_title('Predicted Profit vs. Population Size')
plt.show()

由于梯度方程式函数也在每个训练迭代中输出一个代价的向量，所以我们也可以绘制。 请注意，代价总是降低 - 这是凸优化问题的一个例子。

fig, ax = plt.subplots(figsize=(12,8))
ax.plot(np.arange(iters), cost, 'r')
ax.set_xlabel('Iterations')
ax.set_ylabel('Cost')
ax.set_title('Error vs. Training Epoch')
plt.show()

我们也可以使用scikit-learn的线性回归函数，而不是从头开始实现这些算法。 我们将scikit-learn的线性回归算法应用于第1部分的数据，并看看它的表现。

# set X (training data) and y (target variable)
path =  'linear_regression.txt'
data = pd.read_csv(path, header=None, names=['Population', 'Profit'])

cols = data.shape[1]
X = data.iloc[:,:cols-1]#X是所有行，去掉最后一列
y = data.iloc[:,cols-1:cols]#X是所有行，最后一列

from sklearn import linear_model
model = linear_model.LinearRegression()
model.fit(X, y)
# 输出回归系数
print("回归系数：", model.coef_)
print("截距：", model.intercept_)

回归系数： [[1.19303364]]
截距： [-3.89578088]

print(X.shape)
print(y.shape)

(97, 2)
(97, 1)

scikit-learn model的预测表现

from sklearn.metrics import mean_squared_error, r2_score
y_pred = model.predict(X)
# 评估模型性能
mse = mean_squared_error(y, y_pred)
r2 = r2_score(y, y_pred)
print("均方误差(MSE)：", mse)
print("R平方值：", r2)

# 绘制训练集和测试集的回归线
plt.figure(figsize=(10, 6))
plt.scatter(X, y, color='blue', label='Traning Data')
plt.plot(X, model.predict(X), color='red', linewidth=2, label='Prediction')
plt.xlabel('Population')
plt.ylabel('Profit')
plt.title('Predicted Profit vs. Population Size')
plt.legend()
plt.grid(True)
plt.show()

均方误差(MSE)： 8.953942751950358
R平方值： 0.7020315537841397