机器学习：简单的随机梯度下降(SGD)求解回归问题和用逻辑回归(LR)解决分类问题

20210405作业：

1. 回归问题: 用随机梯度下降法实现，数据用data.csv。
2. 分类问题: 用梯度下降实现逻辑回归，可以用批量梯度也可以用随机梯度实现。数据采用西瓜数据3.0α.csv。

1. SGD (Stochastic gradient descent)

# 导包
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# 导文件
data = pd.DataFrame(pd.read_csv('数据集/data.csv'))
data

sgd:

1. np.random.randn()取随机数组初始化w（符合正态分布），代替[0,0]，在后续学习过程中更加合理 吴恩达机器学习——随机初始化（random initialization）
2. W记录前一次j循环的结果，_W为当前j循环内时时更新的wj的储存数组，待所有wj全部更新完成后，用_W更新W。当前正在梯度更新的wj记入_W
3. 每次进入i循环，初始化cost代价函数，最小二乘法计算代价 最小二乘法本质和原理（least sqaure method）。当前代价与上一代价pre_cost相减，直到逼近其最小值，即小于某一域值时，跳出i循环，返回W
4. 和BGD很像，SGD交换i循环内，j和k循环的先后（当然，其后range也要同步替换）。
   BGD通过遍历所有点，依次更新wj；而SGD遍历一点，更新全部wj，再遍历下一点…所以，SGD可能在未遍历完所有点的情况下，找到符合设定域值的一组W，速度更快，但准确性有所下降
5. 梯度下降就是一种根据代价函数的导数斜率去不断更改系数W，从而达到使导数接近于0，也就是代价函数达到极值点的一种算法

def sgd(x,y,alpha,n):
    """
    sgd函数，求w
    :param x: 数据集的x属性list
    :param y: 数据集的y标签list
    :param alpha: 梯度下降学习率
    :param n: 循环次数
    :return: w随机梯度下降所求系数list
    """
    # w
    w = np.random.randn(2)
    _w = np.empty(2, dtype=float)
    _w = w
    # 数据集样本数m，属性数d
    m = x.shape[0]
    d = w.shape[0]
    # 代价
    pre_cost = 0

    for i in range(n):
        for k in range(m):
            for j in range(d):
                _w[j] = _w[j]-alpha*(w[0]*x[k][0]+w[1]*x[k][1]-y[k])*x[k][j]
        for j in range(d):
            w[j] = _w[j]
        cost = 0
        for k in range(m):
            cost += (w[0]*x[k][0]+w[1]*x[k][1]-y[k])**2
        cost = cost/m
        if abs(cost-pre_cost)<0.001:
            break;
        pre_cost=cost
    return w

注：

1. pandas的.values以array形式返回制定column的取值

# 构造x列表
x_list = []
x_arr = data['x'].values
for i in range(x_arr.shape[0]):
    temp = [1,x_arr[i]]
    x_list.append(temp)
x_list = np.array(x_list)
x_list

# 构造y列表
y_list = np.array(data['y'].values)
y_list

# 计算w
w = sgd(x_list,y_list,alpha=0.00000005,n=10000)
w

注：

1. np.linspace(start, stop, num, endpoint=True, retstep=False, dtype=None)：
   在指定的start和stop，返回固定间隔num的数据
2. plt.scatter(X[], Y[])：“撒点”

# 绘图
x = np.linspace(400,900,100)
y = w[1]*x + w[0]
plt.title('Graph of SGD')
plt.plot(x,y,color='green')
plt.scatter(data['x'].values, data['y'].values)
plt.grid()
plt.xlabel('x')
plt.ylabel('y')
plt.show()

---

2. LR (Logistic Regression)

# 导数据，并简单查看一下，发现没有表头，所以加上表头x1，x2和y
xg_data = pd.DataFrame(pd.read_csv('数据集/西瓜数据3.0α.csv', header=None, names=['x1', 'x2', 'y']))
xg_data.head(5)

sigmoid：

1. S形曲线，g(z)=1/(1+exp(-z))，表示落在某值的概率
2. z是np.mat类型的，那么需要用np.exp而不是math.exp()
3. Sigmoid曲线的’中点’，也就是g(z)=0.5那点，表示概率是50%，将正反例分隔开，
   令z=0，即WX+b=0，可求得决策边界（Decision Boundary），
   plot_x2 = -w[0]/w[2]-(w[1]/w[2])plot_x1

Sigmoid函数概念
机器学习中 Sigmoid 函数的物理含义

LR：

1. 这里求W采用BGD矩阵法
2. 这里的代价函数cost本质上是通过极大似然估计得出的，在上文课件里有具体步骤

极大似然估计原理(Maximum Likelihood Estimate)

def sigmoid(z):
    return 1/(1+np.exp(-z))

def LR(x, y, alpha, n):
    """
    计算逻辑回归的w
    :param x: 样本属性
    :param y: 样本标签
    :param alpha: 学习率
    :param n: 循环次数
    :return: sigmoid曲线的变量w
    """
    # 样本数m
    m = x.shape[0]
    w = np.mat(np.random.randn(x.shape[1],1))
    pre_cost = 0
    for i in range(n):
        f = sigmoid(x*w)
        w = w + alpha*x.T*(y-f)
        cost = (y.T*np.log(f)+(1-y).T*np.log(1-f))/m
        if abs(pre_cost-cost)<0.0001:
            break
        pre_cost = cost
    return w

注：

1. zip将两list合并成tuple列表，从而构造X矩阵。
   这里因为zip在python新版本中返回的是地址，而要显示得用list()转换成列表。
   当然也可以用[:,:]来切片，都可以啦
2. np.insert(arr,obj,values,axis=None)：arr原matrix，
   obj插入位置，values插入新矩阵值，axis若不设置，则会降维
3. np.reshape()中，后参数[x,y]指定重塑成x行y列的ndarray，
   若x或y为-1，则表明在满足另一维度的情况下，自动调整当前维度元素个数
4. (matrix).getA()返回matrix的数组，因为matplotlib作图时，需要array而不是matrix型数据
5. 作图时，用到了布尔去筛选所需样本点

# 获取样本属性X的矩阵
x = list(zip(list(xg_data['x1']), list(xg_data['x2'])))
x = np.mat(x)
one_mat = np.ones(x.shape[0])
x = np.insert(x, 0,values=one_mat, axis=1)
x

# 获取样本标签Y的矩阵
y = np.mat(np.reshape(np.asarray(xg_data['y']), [-1,1]))
y

# 求w
w = LR(x, y, alpha=0.5, n=100000)
w = w.getA()
w

# 结果展示
bl_mat = xg_data['y']==1
plt.scatter(xg_data['x1'][bl_mat], xg_data['x2'][bl_mat], c='green', label='Positive')
plt.scatter(xg_data['x1'][bl_mat!=True], xg_data['x2'][bl_mat!=True], c='red', label='Negative')
plt.xlabel('x1')
plt.ylabel('x2')
plot_x1 = np.linspace(0.1, 0.8, 100)
plot_x2 = -w[0][0]/w[2][0]-(w[1][0]/w[2][0])*plot_x1
plt.plot(plot_x1, plot_x2, color='blue', label='Decision Boundary')
plt.grid()
plt.legend()
plt.show()