随机梯度下降底层实现

import numpy as np
from numpy import *
import matplotlib.pyplot as plt
from sklearn import linear_model

'''
随机梯度下降
'''

# 显示中文
from pylab import *  # 显示中文

mpl.rcParams['font.sans-serif'] = ['SimHei']  # 显示中文

# 画图中显示负号
import matplotlib

matplotlib.rcParams['axes.unicode_minus'] = False

# SGD：每一代(epoch/iteration)使所有样本都参加梯度下降1次：每次梯度下降以一个样本计算梯度，并记录和显示每个样本的代价
# 终止条件：epoch/iteration达到一定条件数目，或者两次theta差值小于某一值

# 数据预处理
# 数据加载
data = np.loadtxt('ex1data1.txt', delimiter=',');
print(data)
# order=np.random.permutation(data.shape[0])
# data=data[order]
# 提取数据X,y
X = data[:, :-1];
print(X)
y = data[:, -1];
print(y)

# 数据的初始化/标准化
m = X.shape[0]  # 样本个数m
X = np.c_[np.ones(m), X];
print(X)
y = np.c_[y];
print(y)


# 定义模型model
def model(X, theta):
    return np.dot(X, theta)


# 定义代价函数costFunction
def costFunction(h, y):
    m = y.shape[0]
    J = 1.0 / (2 * m) * np.dot((h - y).T, (h - y))  # J.shape=() 当h.shape=(1,1),y.shape=(1,1)
    # 或 J = 1.0 / (2.0 * m) * np.sum((h - y) * (h - y))
    return J


# 计算精度的函数
def score(X, y, theta):
    h = model(X, theta)
    y_mean = np.mean(y)
    u = np.sum((h - y) ** 2)  # 计算预测值与真实值的方差
    v = np.sum((y - y_mean) ** 2)  # 计算真实值方差
    score = 1 - u / v
    return score


# 定义梯度下降算法gradDesc
def gradDesc(X, y, alpha=0.005, epoch=150, epsilon=1e-330):
    m, n = X.shape  # 样本个数m，列数n
    theta = np.zeros((n, 1))  # 初始化theta为0
    J_history = []  # 初始化历史代价记录

    theta_old = theta  # 初始化theta前值
    count = 0  # 记录迭代步数

    # 随机梯度下降算法
    for i in range(epoch):

        # 每次全部样本训练完后，重新洗牌  随机梯度下降，随机选取样本
        order = np.random.permutation(m)  # 随机梯度下降，随机选取样本  每一次epoch的不一样
        X, y = X[order], y[order]

        # 每次选取一个样本进行梯度下降
        for j in range(m):
            # print('=========================================')
            # print(X[j])
            # print(X[j:j + 1])
            # print('=========================================')

            h = model(X[j], theta)  # 预测值
            J_history.append(costFunction(h, y[j]))  # 记录代价

            # 计算梯度，并更新theta    选一个原来1/m  m现在是1
            delta_theta = 1.0 * np.dot(X[j:j + 1].T, (h - y[j:j + 1]))  # 计算梯度:x[j:j+1]和[j]不同:前者返回二维，后者一维
            # delta_theta=1.0/m*(X[j]*(h-y[j])).reshape(-1,1)  #计算梯度
            theta -= alpha * delta_theta  # 更新theta
            count += 1  # 记录迭代次数

            # 判断是否终止  norm范数 默认是2 2范数
            if np.linalg.norm(theta - theta_old) < epsilon:  # 若满足两次theta差小于epsilon
                break
            else:
                theta_old = theta

    return theta, J_history, count


# 调用梯度下降算法
theta, J_history, count = gradDesc(X, y)
print(theta)
print('count=', count)
# 画代价曲线新
plt.figure('代价曲线')
plt.title('代价曲线')
plt.xlabel('迭代次数')
plt.ylabel('代价')
print(J_history);
print(np.array(J_history))
# squeeze 函数：从数组的形状中删除单维度条目，即把shape中为1的维度去掉
plt.plot(np.squeeze(J_history))  # np.array(J_history).shape(150,1,1)
plt.show()

# 计算精度并输出
score = score(X, y, theta)
print('精度=', score * 100, '%')

# 画出所有的样本及回归直线
plt.figure('样本及回归直线')
plt.title('样本及回归直线：精度=' + str(round(score, 3)))
plt.scatter(X[:, 1], y[:, 0], c='r', marker='x')
h = model(X, theta)
plt.plot(X[:, 1], h[:, 0])
plt.xlabel('X1')
plt.ylabel('X2')
plt.show()