多元统计分析 矩阵代数

Question

多元统计分析矩阵代数题目求解，讲解下面图片里的题目，写一下解题过程，非常感谢

Answer 1

#!/usr/bin/env python
# -*- coding: utf-8 -*-
# @Time : 2022/4/5 20:52
# @Author : cc
# @File : multi-variable-gd.py
# @Software: PyCharm
# 房价预测。 ex1data2.txt：面积、卧室数、房价
import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
 
# 代价函数，返回该模型w，b参数下的代价
def cost_function(w_hat_matrix, X_matrix, y_matrix) -> float:
    """ 此方法不管几维都固定不变
    :param w_hat_matrix: 二维matrix，1 * 3
    :param X_matrix: 二维matrix，47 * 3
    :param y_matrix: 二维matrix，47 * 1
    :return:  代价
    """
    m = len(X_matrix)  # 样本数
    return np.sum(np.power(X_matrix * w_hat_matrix.T - y_matrix, 2)) / (2 * m)
 
# 梯度下降函数
def gradientDescent(alpha, iters, w_hat_matrix, X_matrix, y_matrix) -> tuple:
    """
    :param alpha: 梯度下降学习率/步长
    :param iters:   梯度下降次数
    :param w_hat_matrix:   二维matrix，1 * 3；
                    一般初始设置为[[0, 0, 0]]，初值对梯度下降收敛速度影响大
    :param X_matrix: 二维matrix，47 * 3。就是公式推导里记的那个超大矩阵X.存放扩展了1的所有样本
    :param y_matrix: 二维matrix，47 * 1。存放真实标记
    :return:
    """
    parameters = int(w_hat_matrix.shape[1])  # 3个参数
    m = len(X_matrix)  # 样本数47
    cur_w_mat_matrix = np.matrix(np.zeros(w_hat_matrix.shape))  # 暂存每次迭代得到的w_hat1*3
    cost = np.zeros(iters)  # 记录每次迭代后的新的代价cost
    # 迭代iters次
    for i in range(iters):
        # 所有样本x_i预测输出 和真实标记 y_i的误差。第i行记录了样本x_i的误差
        error_matrix = X_matrix * w_hat_matrix.T - y_matrix  # 97*1
        # 梯度下降公式：得到新的w_hat：j指向每一列，进行更新（w_1, w_2,..,w_j,...w_d, b）
        for j in range(parameters):
            # np.multiply：m*n和m*n的相同下标元素元素相乘，结果还是m*n的矩阵。可以看成多个数乘
            s = np.sum(np.multiply(error_matrix, X_matrix[:, j]))  # 公式：error矩阵和x的第j列相乘
            cur_w_mat_matrix[0, j] = cur_w_mat_matrix[0, j] - alpha * s / m
        w_hat_matrix = cur_w_mat_matrix  # 每次下降更新w_hat
        cost[i] = cost_function(w_hat_matrix, X_matrix, y_matrix)  # 更新代价，用于观察梯度下降的代价变化曲线
    return w_hat_matrix, cost
 
# 画结果图
def plot_res(transform_w_hat, data_origin, mins, maxs):
    """
    :param transform_w_hat: 符合缩放前原始样例的 b,w1,w2,...wm
            二维array：[[88307.21151185   138.22534685 -7709.05876589]]
    :param data_origin: 原始数据每列的值：也是二维array形式
            属性x1列的所有值array：data_origin[:, 0]
            属性x2...            data_origin[:, 1]
    :param mins:    原始数据每列最小值
                    array: [   852      1 169900]
    :param maxs:    原始数据每列最大值
                    array: [  4478      5 699900]
    :return:  画出拟合模型图和散点图
    """
    # 建立三维模型
    fig = plt.figure() # Create a new figure, or activate an existing figure.
    ax = Axes3D(fig, auto_add_to_figure=False) # 三维坐标轴
    fig.add_axes(ax) # 给图形fig添加坐标轴 an Axes to the figure.
    # 设置三维图角度
    ax.view_init(elev=25, azim=125) # 10 80观察更好
    # 设置三根轴的名称
    ax.set_xlabel('Size')
    ax.set_ylabel('Bedrooms')
    ax.set_zlabel('Prices')
    # 设置x1 x2轴范围
    x1 = np.arange(mins[0], maxs[0] + 1, 1) # x1轴的范围：步长为1
    x2 = np.arange(mins[1], maxs[1] + 1, 1) # x2轴的范围：步长为1
    x1, x2 = np.meshgrid(x1, x2) # 生成网格点坐标矩阵，这句话必须有
 
    # 画线性回归模型：平面图
    b, w1, w2 = transform_w_hat[0, 0], transform_w_hat[0, 1], transform_w_hat[0, 2]# 获取系数
    f = b + w1 * x1 + w2 * x2     # 模型：映射关系
    ax.plot_surface(x1, x2, f, rstride=1, cstride=1, color='red') # 创建平面图（模型）
    # 创画样例散点图
    ax.scatter(data_origin[:, 0], data_origin[:, 1], data_origin[:, 2])
    plt.show()
 
# 每次梯度下降的代价变化图
def plot_cost(cost, iters: int):
    """
    :param cost: 一维array，第i个元素存放第i次梯度下降时的代价
    :param iters: 迭代次数，固定1000次
    :return:
    """
    # 二维坐标轴直接Plt：设置坐标轴名称和标题
    plt.xlabel("iterations")
    plt.ylabel("Cost")
    plt.title("Error vs Traning Epoch")
    # 画直线，x变化范围为0~迭代次数，y为每次的代价
    plt.plot(range(iters), cost, color='red')
    plt.show()
 
# 把数据经过特征缩放（均值标准化）的w_hat变成 符合原始数据的w_hat_transform
def w_hat_transform(arr_w_hat_T, means_T, stds_T):
    """ 同型array相乘相除：对应位置元素相乘，返回矩阵仍然是原型
    :param arr_w_hat:  3*1 的二维array
                array([[-1.11069546e-16],
                    [ 8.78503652e-01],
                    [-4.69166570e-02]])
    :param means_T: 3*1 的二维array
    :param stds_T: 3*1 的二维array
    :return:  咱也不知道为啥这么缩放，抄就完事了。。。。
        标准化的公式 ： data = (data - data.mean()) / data.std()
        转化：
            1. temp = y的均值  * w / y的标准差
            2. 转化的b = (b - sum(temp)) * y的标准差  + y的均值
            3. 转化的w = w * y的标准差 / x的标准差
        最后把w_hat_T恢复成 1 * 3返回
    """
    # data = (data - data.mean()) / data.std()
    temp = means_T[:-1] * arr_w_hat_T[1:] / stds_T[:-1]
    arr_w_hat_T[0] = (arr_w_hat_T[0] - np.sum(temp)) * stds_T[-1] + means_T[-1]
    arr_w_hat_T[1:] = arr_w_hat_T[1:] * stds_T[-1] / stds_T[:-1]
    return arr_w_hat_T.reshape(1, -1)
 
if __name__ == '__main__':
    """ pandas得到dataframe类型的数据data    
            Size    Bedrooms    Price  3列
    """
    path = 'ex1data2.txt'
    data = pd.read_csv(path, names=["Sizes", "Bedrooms", "Prices"])
    # 获取原始数据data的一些描述
    data_origin = data.values  # 二维array，每个子array存放样例的数据
    means = data.mean().values # 一维array，第i个元素：第i列的 均值
    stds = data.std().values  # 一维array，第i个元素：第i列的 标准差
    mins = data.min().values  # 一维array，第i个元素：第i列的 最小值
    maxs = data.max().values # 一维array，第i个元素：第i列的 最大值
    '''data特征缩放：均值标准化，让不同特征值差异不要太大，否则梯度下降收敛会很慢'''
    data = (data - data.mean()) / data.std()
 
    '''添加列，用于获得x_hat组成的X：详情见公式推导，为了计算省略b
            Ones    Size    Bedrooms    Price  4列
    '''
    data.insert(0, 'Ones', 1)  # 在第0列，插入一列属性值全为1的列，列名Ones
    cols = data.shape[1]  # 列数4
    """对于dataframe对象使用iloc
    [p1, p2]
    p1:表示取哪些行：  
        a:b：提取行a~b-1
        [c]：提取行c
    p2:表示取哪些列
        a:b：提取列a~b-1
        [c]：提取列c
    输出形式都是：dataframe的形式
    """
    data_x_hat = data.iloc[:, 0: cols - 1]  # 取dataframe的x：前3列
    data_y = data.iloc[:, [cols - 1]]  # 取dataframe 的y
    # 获取matrix类型的所有样本X和真实标记y_matrix：有的计算只能二维matrix做
    X = np.matrix(data_x_hat) # 47 * 3
    y_matrix = np.matrix(data_y) # 47 * 1
    w_hat_matrix = np.matrix([0, 0, 0])  # 1*3 初始w一般设置为0
    # 看看获取的matrix是几行几列： (47, 3) (47, 1) (1, 3)
    # print(X.shape, y_matrix.shape, w_hat_matrix.shape)
    # 设置步长和下降次数
    alpha = 0.01 # 常考率 0.01， 0.03， 0.1， 0.3， 1， 3， 10
    iters = 1000  # 迭代次数
    # 梯度下降iters次后，我们获得的res_w_hat能使得代价函数获得【局部最小值】
    res_w_hat, cost = gradientDescent(alpha, iters, w_hat_matrix, X, y_matrix)
    print("res_w_hat：", res_w_hat)
    # 用能使得代价函数获得【局部最小值】的res_w_hat，带入计算局部最小代价
    res_cost = cost_function(res_w_hat, X, y_matrix)
    print("res_cost：", res_cost)
 
    """之前获得的res_w_hat，是数据缩放后得到的res_w_hat
    最终我们要把res_w_hat恢复成与原始数据对应的transform_w_hat
    data = (data - data.mean()) / data.std()
    """
    """reshape(a, b)
    如：np.matrix w = [[-1.11069546e-16  8.78503652e-01 -4.69166570e-02]]
    a/b = -1表示不关心行数/列数
    w.reshape(-1,1) 表示把1行3列的matrix变成1列的matrix，形式为
    matrix([[-1.11069546e-16],
        [ 8.78503652e-01],
        [-4.69166570e-02]])
    pandas变成excel就是只有一列    
    """
    # 咋说呢，就是做了个转置
    res_w_hat_reshape = np.array(res_w_hat.reshape(-1, 1)) # 把1*3的二维matrix变成 3*1 的二维array
    means_reshape = means.reshape(-1, 1)  # 3 * 1 的二维array
    stds_reshape = stds.reshape(-1, 1) # 3*1 的二维array
    transform_w_hat = w_hat_transform(res_w_hat_reshape, means_reshape, stds_reshape)
    print("transform_w_hat:", transform_w_hat)
    plot_res(transform_w_hat, data_origin, mins, maxs) # 模型
    plot_cost(cost, iters) # 代价函数曲线