机器学习练习题3.3编程实现对率回归

水热釜炼丹师

于 2023-03-03 11:51:15 发布

阅读量286

点赞数 1

分类专栏：机器学习文章标签： python 机器学习 Powered by 金山文档

本文链接：https://blog.csdn.net/pppgo/article/details/129317282

版权

机器学习专栏收录该内容

1 篇文章 0 订阅

订阅专栏

参考代码

本文链接：https://blog.csdn.net/huzimu_/article/details/123760239

########################################################################我做了什么？

对代码进行部分修改并注解，将查询到的代码解释以及链接附上，从0学习记录

########################################################################

代码、结果、csv

import numpy as np
import pandas as pd

import matplotlib.pyplot as plt




from sklearn import linear_model


# Sigmoid 函数 p58 形似S的连续函数
def sigmoid(s):
    s = 1 / (1 + np.exp(-s))  #3.17
    return s
# 3.27
def J_cost(X, y, beta):
    """

    :param X: array , shape(num_sample, num_features)
    :param y: array, shape(num_sample, 1)
    :param beta: array, shape(num_features+1, 1)
    :return: 公式3.27的结果，即对数似然相反数的值
    """
    X_hat = np.c_[X, np.ones((X.shape[0], 1))]  # 公式3.9上方矩阵
    # mp.c_ 用于连接两个矩阵 https://blog.csdn.net/qq_33728095/article/details/102512600
    # ones函数用于创建指定长度或形状的全1数组，np.ones(3,1) 创建三行一列的全1数组 https://blog.csdn.net/weixin_44884379/article/details/112171669
    # shape[0]读取矩阵第一维度的长度，即数组的行数 https://blog.csdn.net/wzk4869/article/details/126022909
    y = y.reshape(-1, 1)  # 改成m行n列的矩阵格式，-1表示自动计算行或列数  https://blog.csdn.net/a8039974/article/details/119925040?ops_request_misc=%257B%2522request%255Fid%2522%253A%2522167774109816782427446849%2522%252C%2522scm%2522%253A%252220140713.130102334..%2522%257D&request_id=167774109816782427446849&biz_id=0&utm_medium=distribute.pc_search_result.none-task-blog-2~all~top_positive~default-1-119925040-null-null.142^v73^insert_down4,201^v4^add_ask,239^v2^insert_chatgpt&utm_term=reshape&spm=1018.2226.3001.4187
    beta = beta.reshape(-1, 1)
    lbeta = -y*np.dot(X_hat, beta) + np.log(1 + np.exp(np.dot(X_hat, beta)))
    # np.dot响亮的点积和矩阵乘法 https://blog.csdn.net/Liang_xj/article/details/85003467?ops_request_misc=%257B%2522request%255Fid%2522%253A%2522167773576816800192253961%2522%252C%2522scm%2522%253A%252220140713.130102334..%2522%257D&request_id=167773576816800192253961&biz_id=0&utm_medium=distribute.pc_search_result.none-task-blog-2~all~top_positive~default-1-85003467-null-null.142^v73^insert_down4,201^v4^add_ask,239^v2^insert_chatgpt&utm_term=np.dot&spm=1018.2226.3001.4187
    return lbeta.sum()

# 公式3.30 beta的一阶导数


def gradient(X, y, beta):
    """

    :param X:
    :param y:
    :param beta:
    :return:
    """
    X_hat = np.c_[X, np.ones((X.shape[0], 1))]
    y = y.reshape(-1, 1)
    beta = beta.reshape(-1, 1)
    p1 = np.exp(np.dot(X_hat, beta)) / (1 + np.exp(np.dot(X_hat, beta)))
    grad = (-X_hat*(y-p1)).sum(0)
    # sum(0) 每一列求和 https://blog.csdn.net/maple05/article/details/107943144?spm=1001.2101.3001.6650.1&utm_medium=distribute.pc_relevant.none-task-blog-2%7Edefault%7ECTRLIST%7ERate-1-107943144-blog-118693790.pc_relevant_3mothn_strategy_recovery&depth_1-utm_source=distribute.pc_relevant.none-task-blog-2%7Edefault%7ECTRLIST%7ERate-1-107943144-blog-118693790.pc_relevant_3mothn_strategy_recovery&utm_relevant_index=2
    return grad.reshape(-1, 1)

# 梯度下降，更新beta


def update_parameters_gradDesc(X, y, beta, learning_rate, num_iterations, print_cost):
    """

    :param X:
    :param y:
    :param beta:
    :param learning_rate:
    :param num_iterations:
    :param print_cost:
    :return:
    """
    for i in range(num_iterations):
        grad = gradient(X, y, beta)
        beta -= learning_rate * grad

        if (i % 100 == 0) & print_cost:
            print('{}th iteration, cost is {}'.format(i, J_cost(X, y, beta)))
    return beta

# 牛顿法，海森因矩阵


def hessian(X, y, beta):
    """
    根据公式3.31求lbeta二阶导
    :param X:
    :param y:
    :param beta:
    :return:
    """
    X_hat = np.c_[X, np.ones((X.shape[0], 1))]
    y = y.reshape(-1, 1)
    beta = beta.reshape(-1, 1)
    p1 = np.exp(np.dot(X_hat, beta)) / (1 + np.exp(np.dot(X_hat, beta)))
    m, n = X.shape
    P = np.eye(m) * p1 * (1 - p1)  # np.eye(m)为单位阵，P为方阵
    assert P.shape[0] == P.shape[1]  # assert断言函数，为真继续执行，为假报错 https://blog.csdn.net/weixin_61561736/article/details/124886522?ops_request_misc=%257B%2522request%255Fid%2522%253A%2522167774221816782425188373%2522%252C%2522scm%2522%253A%252220140713.130102334..%2522%257D&request_id=167774221816782425188373&biz_id=0&utm_medium=distribute.pc_search_result.none-task-blog-2~all~top_positive~default-1-124886522-null-null.142^v73^insert_down4,201^v4^add_ask,239^v2^insert_chatgpt&utm_term=assert&spm=1018.2226.3001.4187
    hess = np.dot(np.dot(X_hat.T, P), X_hat)  # .T 转置
    return hess

# 牛顿法更新参数


def update_parameter_Newton(X, y, beta, num_iterations, print_cost):
    """
    根据公式3.29更新beta
    :param X:
    :param y:
    :param beta:
    :param num_iterations:
    :param print_cost:
    :return:
    """
    for i in range(num_iterations):
        grad = gradient(X, y, beta)
        hess = hessian(X, y, beta)
        beta -= np.dot(np.linalg.inv(hess), grad)  # np.linalg.inv矩阵求逆

        if (i % 100 == 0) & print_cost:
            print('{}th iteration, cost is {}'.format(i, J_cost(X, y, beta)))
    return beta

# 初始化beta


def initialize_beta(n):
    beta = np.random.randn(n + 1, 1) * 0.5 + 1
    # randn产生的样本符合正态分布 https://blog.csdn.net/u012149181/article/details/78913167?ops_request_misc=%257B%2522request%255Fid%2522%253A%2522167774570616800222825721%2522%252C%2522scm%2522%253A%252220140713.130102334..%2522%257D&request_id=167774570616800222825721&biz_id=0&utm_medium=distribute.pc_search_result.none-task-blog-2~all~top_positive~default-1-78913167-null-null.142^v73^insert_down4,201^v4^add_ask,239^v2^insert_chatgpt&utm_term=np.random.randn&spm=1018.2226.3001.4187
    return beta

# 对数几率模型


def logistic_beta(X, y, num_interation=100, learning_rate=1.2, print_cost=False, method='gradDesc'):
    """

    :param X: array shape(num_sample, num_feature)
    :param y:
    :param num_interation: 迭代次数
    :param learning_rate: 学习率
    :param print_cost: 是否打印对数似然相反数
    :param method: 所用方法
    :return:
    """
    m, n = X.shape
    beta = initialize_beta(n)

    if method == 'gradDesc':
        return update_parameters_gradDesc(X, y, beta, learning_rate, num_interation, print_cost)
    elif method == 'Newton':
        return update_parameter_Newton(X, y, beta, num_interation, print_cost)
    else:
        return ValueError('Unknown error of %s' % method)

#预测


def predict(X, beta):
    X_hat = np.c_[X, np.ones((X.shape[0], 1))]
    beta.reshape(-1, 1)
    p1 = np.exp(np.dot(X_hat, beta)) / (1 + np.exp(np.dot(X_hat, beta)))
    p1[p1 > 0.5] = 1
    p1[p1 <= 0.5] = 0
    temp = np.c_[X, p1]
    print(temp)
    return p1

if __name__ == '__main__':
    data_path = r'G:\pythonProject\machinelearning\watermelon3.csv'
    # 加载数据
    data = pd.read_csv(data_path, encoding='ISO-8859-1').values
    # print(data)
    is_good = data[:, 2] == 'yes'  # : 表示全部行的数据 9-第九列
    is_bad = data[:, 2] == 'no'
    # 7，8列数据，密度、含糖量
    X = data[:, 0:2].astype(float)
    y = data[:, 2]
    y[y == 'yes'] = 1
    y[y == 'no'] = 0
    y = y.astype(int)

    # 绘图：密度为行，含糖量为列，好瓜为黑圈，坏瓜为红x
    plt.scatter(data[:, 0][is_good], data[:, 1][is_good], c='k', marker='o')
    plt.scatter(data[:, 0][is_bad], data[:, 1][is_bad], c='r', marker='x')

    plt.xlabel('density')
    plt.ylabel('sugar')

    # 可视化模型结果gradDesc
    beta = logistic_beta(X, y, num_interation=1000, learning_rate=2, print_cost=True, method='gradDesc')
    w1, w2, intercept = beta
    #预测结果
    predict(X, beta)
    x1 = np.linspace(0, 1)
    # 创建数值序列 https://blog.csdn.net/neweastsun/article/details/99676029?ops_request_misc=%257B%2522request%255Fid%2522%253A%2522167774351316800215063759%2522%252C%2522scm%2522%253A%252220140713.130102334..%2522%257D&request_id=167774351316800215063759&biz_id=0&utm_medium=distribute.pc_search_result.none-task-blog-2~all~top_positive~default-1-99676029-null-null.142^v73^insert_down4,201^v4^add_ask,239^v2^insert_chatgpt&utm_term=np.linspace&spm=1018.2226.3001.4187
    y1 = -(w1 * x1 + intercept) / w2
    plt.plot(x1, y1, label=r'my_logistic_gradDesc')

    # 可视化模型结果gradDesc
    beta = logistic_beta(X, y, num_interation=1000, learning_rate=2, print_cost=True, method='Newton')
    w1, w2, intercept = beta
    #预测结果
    predict(X, beta)
    x1 = np.linspace(0, 1)
    y1 = -(w1 * x1 + intercept) / w2
    plt.plot(x1, y1, label=r'my_logistic_Newton')

    # # 可视化sklearn Logistic regression 模型结果
    # # 注意sklearn的逻辑回归中，C越大表示正则化程度越低
    # Ir = linear_model.LogisticRegression(solver='lbfgs', c=100)
    # Ir.fit(X,y)
    #
    # Ir_beta = np.c_[Ir.coef_, Ir.intercept_]
    # print('cost of sklearn logistic: {}'.format(J_cost(X, y, beta)))
    #
    # w1_sk, w2_sk = Ir.coef_[0, :]
    #
    # x2 = np.linspace(0, 1)
    # y2 = -(w1_sk * x2 + Ir.intercept_) / w2_sk
    #
    # plt.plot(x2, y2, label=r'sklearn logistic')

    plt.legend(loc='upper right')  # 图例位置
    plt.show()