李宏毅机器学习2020-作业2：判断年收入是否高于50000美元（logistic regression/generative）_通过逻辑回归判别一个人的收入是否超过五万美元-CSDN博客

import numpy as np
import matplotlib.pyplot as plt  # 画图的包= =
# seed里面的参数相同时，生成的随机数也一样
np.random.seed(0)

with open(X_train_fpath) as f:
    next(f)   # 跳过第一行（第一行都是中文解说啦），直接从第二行开始读取数据
    X_train = np.array([line.strip('\n').split(',')[1:] for line in f], dtype=float)  # 列表生成式

with open(Y_train_fpath) as f:
    next(f)
    Y_train = np.array([line.strip('\n').split(',')[1] for line in f], dtype=float)

with open(X_test_fpath) as f:
    next(f)
    X_test = np.array([line.strip('\n').split(',')[1:] for line in f], dtype=float)

3.2 正则化

def _normalize(X, train = True, specified_column = None, X_mean = None, X_std = None):
    '''
       This function normalizes specific columns of X.
       The mean and standard variance of training data will be reused when processing testing data.

       :param X: data to be processed
       :param train: 'True' when processing training data,'False' for testing data
       :param specified_column: indexs of the columns that will be normalized.
                                 if 'none' all column will be normalized.
       :param X_mean: mean value of training data,used when train = 'False'
       :param X_std: standard deviation of training data, used when train = 'False'

       :return:
       X: normalized data
       X_mean:computed mean value of training data
       X_std:computed standard deviation of training data
       '''

    if specified_column == None:
        specified_column = np.arange(X.shape[1])  # arange函数用于创建等差数组

    if train:  # 此时处理的是training data
        X_mean = np.mean(X[:, specified_column], axis=0).reshape(1, -1)
        X_std = np.std(X[:, specified_column], axis=0).reshape(1, -1)

    X[:, specified_column] = (X[:, specified_column] - X_mean) / (X_std + 1e-8)   # 避免除数为0的情况
    return X, X_mean, X_std


X_train, X_mean, X_std = _normalize(X_train, train = True)
X_test, _, _ = _normalize(X_test,train = False,X_mean = X_mean,X_std = X_std)

4. Logistic Regression

4.1 相关公式与推导

对于一组训练集

产生这一组w,b的概率为：

其中

求偏导

最后结果为：

相较于线性回归，异同如下：

4.2 定义需要用到的函数

4.2.1 将训练集分为训练集和验证集

def _train_dev_split(X, Y, dev_ratio = 0.25):
    train_size = int(len(X) * (1 - dev_ratio))
    return X[:train_size], Y[:train_size], X[train_size:], Y[train_size:]
dev_ratio = 0.1
X_train, Y_train, X_dev, Y_dev = _train_dev_split(X_train, Y_train, dev_ratio=dev_ratio)

验证下：

# 求出一些有用的维度数
train_size = X_train.shape[0]
dev_size = X_dev.shape[0]
test_size = X_test.shape[0]
data_dim = X_train.shape[1]

print("Size of training set = {}.".format(train_size))
print("Size of development set = {}.".format(dev_size))
print("Size of test set = {}.".format(test_size))
print("Size of data dim = {}.".format(data_dim))

运行，

4.2.2 行打乱

多维矩阵中，只对第一维（行）做打乱顺序操作。也就是每一行的数据不动，但是行的顺序改变。例：在一维中，np.random.shuffle(randomize) 将列表randomize内元素打乱顺序在二维中，randomize记录X的行下标：randomize = np.arange(len(X))，在经过np.random.shuffle(randomize) 后元素顺序改变；由于randomize与X的下标绑定，randomize内元素顺序改变，那么X的下标也进行同步的改变

def _shuffle(X,Y):
    randomize = np.arange(len(X))
    np.random.shuffle(randomize)
    return X[randomize], Y[randomize]

4.2.3 sigmoid函数

np.clip(a, a_min, a_max, out=None)：

取a数组中的闭区间[a_min, a_max]，数组中小于a_min的数都变成a_min，大于a_max的数都变成a_max

如： a = np.arange(10) np.clip(a, 1, 8) array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9]) 会变成 array([1, 1, 2, 3, 4, 5, 6, 7, 8, 8])

def _sigmoid(z):
    # 为避免溢出，设置了最大最小值，即如果sigmoid函数的最小值比1e-8小，只会输出1e-8；而比1 - (1e-8)大，则只输出1 - (1e-8)
    return np.clip(1 / (1.0 + np.exp(-z)), 1e-8, 1 - (1e-8))

4.2.3 逻辑回归的方程

逻辑回归的方程，输入为x，参数为w，bias是b，注意X与w都是数组，b是一个数

def _f(X,w,b):
    return _sigmoid(np.matmul(X, w) + b)

4.2.3 对预测结果取整数

np.round(数据, decimal=保留的小数位数)
原则：
一般该函数遵循四舍五入原则：np.round(11.5)=12
但是会有特殊情况：当整数部分以0结束时，一律是向下取整：np.round(10.5)=10
比较稳定的浮点数取整：向上的np.ceil(11.5)=12，向下的floor(11.5)=11

def _predict(X, w,b):
    return np.round(_f(X, w, b)).astype(int)

4.2.4 准确率

def _accuracy(Y_pred,Y_label):
    acc = 1 - np.mean(np.abs(Y_pred - Y_label))
    # 如果预测正确，则结果是0，否则结果是1，那么我们求mean平均值的话所得值是1的概率(mean相当于1的个数 / 总个数), 那么我们求0的概率就是1 - it。
    return acc

4.2.5 交叉熵

# 按照公式打出来的
def _cross_entropy_loss(y_pred,Y_label):
    cross_entropy = -np.dot(Y_label, np.log(y_pred)) - np.dot((1 - Y_label), np.log(1 - y_pred))
    return cross_entropy

4.2.6 计算梯度

def _gradient(X,Y_label,w,b):
    y_pred = _f(X, w, b)
    pred_error = Y_lable - y_pred
    # X.T就是X的转置，axis取值为1时代表将每一行的元素相加，实际上返回的是1行510列的数组
    w_grad = -np.sum(pred_error * X.T, 1)
    # 对b求偏微分后的结果，黑板上没有，但因为逻辑回归和线性回归的损失函数相似，可由线性回归对b进行求偏微分得到
    b_grad = -np.sum(pred_error)
    return w_grad, b_grad

4.3 开始训练

# 我们使用小批次梯度下降法来训练。训练资料被分为许多小批次，针对每一个小批次，我们分别计算其梯度以及损失，并根据该批次来更新模型的参数。
# 当一次循环完成，也就是整个训练集的所有小批次都被使用过后，我们将所有的训练资料打散并重新分为新的小批次，进行下一次循环，直至事先设定的循环次数达成位置

# 1.使用0初始化w，b参数
w = np.zeros((data_dim,)) # 默认列数为1
#data_dim = X_train.shape[1]列数
b = np.zeros((1,))  # b就只是一个数字而已啦

max_iters = 20
batch_size = 8
learning_rate = 0.05

# 2.保存每次迭代时的loss和accuracy用于画图
train_loss = []  # 训练集的loss值
dev_loss = []   # 验证集的loss值
train_acc = []   # 训练集的正确率
dev_acc = []    # 验证集的正确率

# 3.迭代训练
# 记录参数更新的次数
step = 1

#Iterative training
for epoch in range(max_iters):
    # 随机的将X，Y的顺序打乱
    X_train, Y_train = _shuffle(X_train, Y_train)

    # Mini_batch training
    for idx in range(int(np.floor(train_size / batch_size))):    # 分别取X和Y中的对应8个数据(每个批次8个数据)
        X = X_train[batch_size * idx: batch_size * (idx + 1)]
        Y = Y_train[batch_size * idx: batch_size * (idx + 1)]

        # 计算梯度
        w_grad, b_grad = _gradient(X, Y, w, b)

        # 更新参数，自适应学习率这次使用的是非常简单的学习率除以更新次数的根
        w = w - learning_rate / np.sqrt(step) * w_grad
        b = b - learning_rate / np.sqrt(step) * b_grad

        step = step + 1

    # 计算training set和devlopment set的loss和准确率
    y_train_pred = _f(X_train, w, b)
    Y_train_pred = np.round(y_train_pred)  # 将数据转换成bool格式的
    train_acc.append(_accuracy(Y_train_pred, Y_train))  # 记录此次迭代的准确率
    train_loss.append(_cross_entropy_loss(y_train_pred, Y_train) / train_size)  #记录此次迭代的loss

    # 对验证集进行同样的操作
    y_dev_pred = _f(X_dev, w, b)
    Y_dev_pred = np.round(y_dev_pred)
    dev_acc.append(_accuracy(Y_dev_pred, Y_dev))
    dev_loss.append(_cross_entropy_loss(y_dev_pred, Y_dev) / dev_size)

# 输出最后依次迭代的结果
print("Training loss : {}.".format(train_loss[-1]))
print("Development loss :{}.".format(dev_loss[-1]))
print("Training acc: {}.".format(train_acc[-1]))
print("Development acc :{}.".format(dev_acc[-1]))
# 将参数保存下来
np.save("Weight.w2.npy", w)

运行，可以看到：

4.4 绘制loss和accuracy图像

# loss curve
plt.plot(train_loss)
plt.plot(dev_loss)
plt.title('loss')
plt.legend(['train', 'dev'])
plt.savefig('loss.png')
plt.show()

# accuracy curve
plt.plot(train_acc)
plt.plot(dev_acc)
plt.title('Accuracy')
plt.legend(['train', 'dev'])
plt.savefig('Acc.png')
plt.show()

4.5 开始测试

# 预测datatest得到预测结果
w = np.load("Weight.w2.npy")  # 加载参数

# Predict testing labels
predictions = _predict(X_test, w, b)
with open(output_fpath.format('logistic'), 'w') as f:
    f.write('id,label\n')
    # # enumerate多用于在for循环中得到计数，利用它可以同时获得索引和值
    for i, label in enumerate(predictions):
        f.write('{},{}\n'.format(i, label))

# # 找到权重中最大的前十项，即关联结果的最紧密的参数
ind = np.argsort(np.abs(w))[::-1]
with open(X_test_fpath) as f:
    content = f.readline().strip('\n').split(',')
features = np.array(content)
for i in ind[0:10]:
    print(features[i], w[i])

5. Porbabilistic generative model 概率生成模型

5.1 计算两个类别的协方差

类别0：年收入不超过50000美元

类别1：年收不超过50000美元

由于两个类别公用一个协方差，所以模型共享的协方差矩阵为所有协方差的加权平均值。

# 第0类为收入高与5w美元的
# 第1类为收入低于5w美元的
X_train_0 = np.array([x for x, y in zip(X_train, Y_train) if y == 0])
X_train_1 = np.array([x for x, y in zip(X_train, Y_train) if y == 1])

mean_0 = np.mean(X_train_0, axis=0)
mean_1 = np.mean(X_train_1, axis=0)

# 计算类内协方差
cov_0 = np.zeros((data_dim, data_dim))
cov_1 = np.zeros((data_dim, data_dim))

# 求第0类的协方差
for x in X_train_0:
    cov_0 += np.dot(np.transpose([x - mean_0]), [x-mean_0]) / X_train_0.shape[0]
# 求第1类的协方差
for x in X_train_1:
    cov_1 += np.dot(np.transpose([x - mean_1]), [x - mean_1]) / X_train_1.shape[0]

print("cov0 : {}".format(cov_0))
print("cov1 : {}".format(cov_1))
#若两类共用一个协方差
# 模型共享的协方差矩阵为所有协方差的加权平均值
cov = (cov_0 * X_train_0.shape[0] + cov_1 * X_train_1.shape[0]) / train_size
print("cov : {}".format(cov))

5.2 求协方差的逆矩阵

np.linalg.svd(a,full_matrices=1,compute_uv=1)
参数：
a是一个形如(M,N)矩阵
full_matrices的取值是为0或者1，默认值为1，这时u的大小为(M,M)，v的大小为(N,N) 。否则u的大小为(M,K)，v的大小为(K,N) ，K=min(M,N)
compute_uv的取值是为0或者1，默认值为1，表示计算u,s,v。为0的时候只计算s。
返回值：
总共有三个返回值u,s,v
u大小为(M,M)，s大小为(M,N)，v大小为(N,N)。
A = u*s*v  qi
其中s是对矩阵a的奇异值分解。s除了对角元素不为0，其他元素都为0，并且对角元素从大到小排列。s中有n个奇异值，一般排在后面的比较接近0，所以仅保留比较大的r个奇异值。

# 计算协方差矩阵的逆
# 协方差矩阵可能是奇异矩阵, 直接使用np.linalg.inv() 可能会产生错误
# 通过SVD矩阵分解，可以快速准确地获得方差矩阵的逆
u, s, v = np.linalg.svd(cov, full_matrices=False)
inv = np.matmul(v.T * 1 / s, u.T)

5.3 参数计算

# 计算w和b
w = np.dot(inv, (mean_0 - mean_1))
b = (-0.5) * np.dot(mean_0, np.dot(inv, mean_0)) + 0.5 * np.dot(mean_1, np.dot(inv, mean_1)) + np.log(float(X_train_0.shape[0]) / X_train_1.shape[0])
# 保存参数
np.save("weight_generative.npy", w)

#加载参数
w = np.load('weight_generative.npy')

Y_trained_pred = 1 - _predict(X_train, w, b)
print('Training accuracy : {}.'.format(_accuracy(Y_trained_pred, Y_train)))

5.4 预测

prediction = _predict(X_test, w, b)
with open(output_fpath.format("generative"),'w') as f:
    f.write("id, label\n")
    for id, label in enumerate(prediction):
        f.write("{}, {}\n".format(id, label))

# 打印出权值最高的10个参数
sort = np.argsort(np.abs(w))[::-1]
with open(X_test_fpath) as f:
    content = f.readline().strip("\n").split(',')
features = np.array(content)
for i in sort[0:10]:
    print(features[i], w[i])

输出的结果为：

结束。