机器学习之线性回归

# 进行数据分析所需库，可以看做是对numpy工具的补充
import pandas as pd
import numpy as np

# 应该把Seaborn视为matplotlib的补充,作图所用工具，在大多数情况下使用seaborn就能做出很具有吸引力的图，而使用matplotlib就能制作具有更多特色的图
import seaborn as sns
import matplotlib.pyplot as plt

# 设置绘画的图标格式和颜色
sns.set(context="notebook", style="whitegrid", palette="dark")
df = pd.read_csv('ex1data1.txt', names=['population', 'profit'])
print(df.head())

sns.lmplot(x='population', y='profit', data=df,
           # size=6,
            fit_reg=False,# 是否画出拟合直线，默认画出
           # scatter_kws={"s": 50}
          )
# 读取数据并赋予列名，此时df共有两列，列名分别为population和profit
plt.show()

def get_x(df):
    # ones:格式化为一个m列的全为1的向量
    ones = pd.DataFrame({'ones': np.ones(len(df))})
    # 将date格式化为ones向量右边加上df矩阵的矩阵
    data = pd.concat([ones, df], axis=1)
    # 返回data矩阵的前两列，返回m * 2的矩阵
    # 这个操作返回 ndarray,不是矩阵
    # Numpy提供了一种存储单一数据类型的多维数组——ndarray
    # 创建一个 ndarray 只需调用 NumPy 的 array 函数即可：
    return data.iloc[:, :-1].as_matrix()


def get_y(df):
    # 返回df矩阵的最后一列，返回m × 1的向量
    return np.array(df.iloc[:, -1])


# def normalize_feature(df):
#     # 特征缩放,对df中的两列数据分别进行特征缩放 ：（x - x平均值）/ x方差
#     return df.apply(lambda column: (column - column.mean()) / column.std())


X = get_x(df)
print(X)
print(X.shape, type(X))  # X为m × 2的矩阵
y = get_y(df)
print(y.shape, type(y))  # y为m × 1的向量
print(y)

theta = np.zeros(X.shape[1])  # X.shape[1]=2,代表特征数n X.shape(0) = 97;X.shape(1) = 2;将theta初始化为2行的零向量（列向量）

# 计算代价函数的值
def lr_cost(theta, X, y):
    #     计算theta固定时，此时的代价函数值
    #     X: R(m*n), m 样本数, n 特征数
    #     y: R(m)
    #     theta : R(n), 线性回归的参数
    m = X.shape[0]  # m为样本数

    # X矩阵和theta矩阵的点积所得矩阵（m × 1） - 矩阵y（m × 1）
    inner = X @ theta - y  # R(m*1)，X @ theta等价于X.dot(theta)

    # square_sum为inner的每个元素平方之和，即二范数
    square_sum = inner.T @ inner

    # 得到此时代价函数的值
    cost = square_sum / (2 * m)

    return cost


# 将代价函数对theta求偏导，再求梯度
def gradient(theta, X, y):
    # m为样本数
    m = X.shape[0]

    # 这里的inner为代价函数对theta求导的结果
    inner = X.T @ (X @ theta - y)  # (m,n).T @ (m, 1) -> (n, 1)，X @ theta等价于X.dot(theta)

    return inner / m


# 梯度下降函数
def batch_gradient_decent(theta, X, y, epoch, alpha=0.01):
    #   拟合线性回归，返回参数和代价
    #     epoch: 批处理的轮数
    #     alpha: theta移动的步长
    #     """

    # 得到初始theta时的代价函数值
    cost_data = [lr_cost(theta, X, y)]

    # 拷贝一份，不和原来的theta混淆
    _theta = theta.copy()

    # 开始迭代epoch次
    for _ in range(epoch):
        # 根据当前梯度更新theta
        _theta = _theta - alpha * gradient(_theta, X, y)
        # 记录每次迭代后代价函数的值
        cost_data.append(lr_cost(_theta, X, y))

    return _theta, cost_data


# 迭代500次求最小代价和所对应的theta
epoch = 500
final_theta, cost_data = batch_gradient_decent(theta, X, y, epoch)
print(final_theta)

# 计算最终的代价函数的值
cost =lr_cost(final_theta, X, y)
print(cost)

# 画出代价函数值变化图
# 可以看到从第二轮代价数据变换很大，接下来平稳了
ax = sns.tsplot(cost_data, time=np.arange(epoch + 1))
ax.set_xlabel('epoch') # 迭代次数
ax.set_ylabel('cost') # 代价函数
plt.show()

b = final_theta[0]  # intercept，Y轴上的截距
m = final_theta[1]  # slope，斜率

# 画出原数据点和线性回归的最终结果
plt.scatter(df.population, df.profit, label="Training data")
plt.plot(df.population, df.population * m + b, label="Prediction")
plt.legend(loc=2)
plt.show()