牛顿方法、指数分布族、广义线性模型和softmax回归—机器学习公开课第四讲

代码实现

牛顿法的代码实现：

# 定义h函数
def h(x1, x2):
    h = theta[0]+theta[1]*x1+theta[2]*x2
    return 1/(1+math.exp(-h))
# 定义迭代次数
cnt = 0
# 定义f函数
def f(x1,x2):
    h = theta[0]+theta[1]*x1+theta[2]*x2
    return 1/((1+math.exp(h))**2)

def newton_method(x1, x2 , y):
    global cnt
    global error0
    while True:
        cnt += 1
        diff = [0]*3
        # 初始化hessian矩阵
        Hessian = [[0, 0, 0], [0, 0, 0], [0, 0, 0]]
        Hessian_mat = np.mat(Hessian,dtype=float)
        # 计算l(theta)导数
        for i in range(m):
            diff[0] += y[i]-h(x1[i], x2[i])
            diff[1] += (y[i]-h(x1[i], x2[i]))*x1[i]
            diff[2] += (y[i]-h(x1[i], x2[i]))*x2[i]
        # 计算hessian矩阵的每一个元素
        for k in range(m):
            Hessian_mat[(0, 0)] += (-f(x1[k], x2[k]))
            Hessian_mat[(0, 1)] += (-f(x1[k], x2[k])) * x1[k]
            Hessian_mat[(0, 2)] += (-f(x1[k], x2[k])) * x2[k]
            Hessian_mat[(1, 0)] += (-f(x1[k], x2[k])) * x1[k]
            Hessian_mat[(1, 1)] += (-f(x1[k], x2[k])) * x1[k] * x1[k]
            Hessian_mat[(1, 2)] += (-f(x1[k], x2[k])) * x1[k] * x2[k]
            Hessian_mat[(2, 0)] += (-f(x1[k], x2[k])) * x2[k]
            Hessian_mat[(2, 1)] += (-f(x1[k], x2[k])) * x2[k] * x1[k]
            Hessian_mat[(2, 2)] += (-f(x1[k], x2[k])) * x2[k] * x2[k]
        Hessian_mat_I = Hessian_mat.I
        theta[0] -= Hessian_mat_I[0]*np.mat(diff).T
        theta[1] -= Hessian_mat_I[1]*np.mat(diff).T
        theta[2] -= Hessian_mat_I[2]*np.mat(diff).T
        error1 = 0
        for i in range(m):
            error1 += (h(x1[i], x2[i]) - y[i]) ** 2 / 2
        print("cnt=%d,Hessian_mat_I=%f,error1=%f" % (cnt, np.linalg.det(Hessian_mat_I), error1))
        if (abs(error1 - error0) < epsilon):
            return (theta,cnt)
        else:
            error0 = error1

softmax回归的代码实现：

import numpy as np
import matplotlib.pyplot as plt

# 迭代次数
iter_num = 3500
# 设置学习率
lr = 0.0001
# 未使用正则化
lamba = 0

# 加载数据集,返回的特征集和标签集都是矩阵
def load_data(filename):
    dataset = []
    labelset = []
    fr = open(filename)
    for row in fr.readlines():
        row = row.strip().split()
        dataset.append([float(row[0]), float(row[1]), float(row[2])])
        labelset.append([int(row[3])])
    return np.mat(dataset), np.mat(labelset)

# 归一化特征集
def normalization(X):
    m = X.shape
    print(m)
    Xmin = np.min(X, axis=0)
    Xmax = np.max(X, axis=0)
    Xmean = np.mean(X, axis=0)
    X_norm = (X-Xmean)/(Xmax-Xmin)
    return X_norm

# sigmoid函数
def sigmoid(x):
    h = 1.0 / (1 + np.exp(-x))
    return h

# 定义一个propab函数
def propab(i, j, x, theta):
    divider_up = np.exp(theta[j] * x[i].T)[0][0, 0]
    divider_down = np.sum(np.exp(theta * x[i].T))
    return divider_up/divider_down

# 定义决策函数
def predict(w,x):
    vector = w*x.T
    y_ = np.argmax(vector)+1
    return y_

def test(w, x, y):
    m1, n1 = x.shape
    cnt = 0
    for i in range(m1):
        y_ = predict(w, x[i])
        if(y_ == y[i, 0]):
            cnt += 1
    accurancy = cnt/m1
    return accurancy

# 定义一个loss函数
def calculate_loss(w, y, x, k):
    loss = 0
    error = 0
    for i in range(m):
        for j in range(k):
            value = propab(i, j, x, w)
            error = int(y[i, 0] == j+1) * np.log(value)
        w_ravel = w.ravel()
        loss += -error/m + lamba * (w_ravel * w_ravel.T)[0, 0] / m
    return loss

# 定义一个softmax回归函数
def softmax_reg(datamat, labelmat, k, Iternum, alpha):
    m, n = datamat.shape
    print("m=", m, " n=", n)
    # 记录正确率和损失值的变化
    accur_list = []
    loss_list = []
    iter_list = []
    # 初始化权重
    w = np.mat(np.ones((k, n)))
    diff = np.mat(np.ones((k, n)))
    iter = 0
    while(iter < Iternum):
        for j in range(k):
            for i in range(m):
                value = propab(i, j, datamat, w)
                diff[j] += -np.mat((int(labelmat[i, 0] == j+1) - value)*datamat[i])/m + lamba * w[j]
            w[j] -= alpha*diff[j]
        iter += 1
        if(iter % 100 == 0):
            print(iter)
            accur = test(w, datamat, labelmat)
            accur_list.append(accur)
            loss = calculate_loss(w, labelmat, datamat, k)
            loss_list.append(loss)
            print("accurancy = ", accur)
            print("loss = ", loss)
            iter_list.append(iter)
    return w, accur_list, loss_list, iter_list

if __name__ == "__main__":
    data_mat, label_mat = load_data("web_yuehui.txt")
    print(data_mat.shape)
    print(label_mat.shape)
    data_norm = normalization(data_mat)
    m, n = data_norm.shape
    offset = np.ones((m, 1))
    trainmat = np.c_[offset, data_norm]
    w_value, accurancy, loss, iter = softmax_reg(trainmat, label_mat, 3, iter_num, lr)
    # 画图
    plt.plot(iter, loss)
    plt.xlabel("iter")
    plt.ylabel("loss")
    plt.savefig("loss_changing.png")
    plt.show()