DataWhale集成学习笔记-task01

DataWhale笔记-task01
学习视频

高等数学复习

1. Hessian矩阵是梯度向量的g(x)对自变量x的雅可比矩阵
   例子：
   $$f(x,y,z)=x^2+y^2+z^2+2x+4y-6z.$$
   梯度向量 -> f(x,y,z)的Jacobian -> g(x,y,z)：
   $$g(x,y,z)={\left(\begin{array}{ccc}2x+2& 2y+4& 2z-6\end{array}\right)}^T.$$
   自变量(x,y,z)的Hessian矩阵 -> g(x,y,z)的Jacobian -> H(x,y,z)
   $$H(x,y,z)=\left(\begin{array}{ccc}2& 0& 0\\ 0& 2& 0\\ 0& 0& 2\end{array}\right).$$
2. 确定极值的重要定律

设n多元实函数 $f\left(x_1,x_2,\cdots ,x_n\right)$ 在点 $M_0\left(a_1,a_2,\dots ,a_n\right)$ 的邻域内有二阶连续偏导，若有:
$${\frac{\partial f}{\partial x_j}|}_{\left(a_1,a_2,\dots ,a_n\right)}=0,j=1,2,\dots ,n$$
并且
$$A=\left[\begin{array}{cccc}\frac{{\partial }^{2}f}{\partial x_1^2}& \frac{{\partial }^{2}f}{\partial x_1\partial x_2}& \cdots & \frac{{\partial }^{2}f}{\partial x_1\partial x_n}\\ \frac{{\partial }^{2}f}{\partial x_2\partial x_1}& \frac{{\partial }^{2}f}{\partial x_2^2}& \cdots & \frac{{\partial }^{2}f}{\partial x_2\partial x_n}\\ ⋮& ⋮& \ddots & ⋮\\ \frac{{\partial }^{2}f}{\partial x_n\partial x_1}& \frac{{\partial }^{2}f}{\partial x_n\partial x_2}& \cdots & \frac{{\partial }^{2}f}{\partial x_n^2}\end{array}\right]$$
则有如下结果:
(1) 当A正定矩阵时, $f\left(x_1,x_2,\cdots ,x_n\right)$ 在 $M_0\left(a_1,a_2,\dots ,a_n\right)$ 处是极小值;
(2) 当A负定矩阵时, $f\left(x_1,x_2,\cdots ,x_n\right)$ 在 $M_0\left(a_1,a_2,\dots ,a_n\right)$ 处是极大值;
(3) 当A不定矩阵时, $M_0\left(a_1,a_2,\dots ,a_n\right)$ 不是极值点。
(4) 当A为半正定矩阵或半负定矩阵时, $M_0\left(a_1,a_2,\dots ,a_n\right)$ 是“可疑"极值点，尚需要利用其他方法来判定。
正定半正定矩阵的定义课见知乎：定义及判断方法

3. 作业：使用梯度下降法求函数 𝑦=𝑐𝑜𝑠(𝑥) 在区间 𝑥∈[0,2𝜋] 的极小值点。

import numpy as np
import matplotlib.pyplot as plt
def fun(x):
    return np.cos(x)
 # 不同的导数算法
def fun_diff1(x):
    return -np.sin(x)
def fun_diff2(f,x,Delta=1e-4):
    return (f(x+Delta)-f(x-Delta))/(2*Delta)
# 梯度下降的方法
def fun_grad(x,grad,max_loop):
    tem=x
    li_x=[x]
    for i in range(max_loop):
        tem=x-fun_diff1(x)*grad
        if x==tem:
            break
        elif (x>2*np.pi) or (x<0):
            print("在0~2pi之间不收敛")
            break
        else:
            x=tem
            li_x.append(x)
    return li_x
init=5
learn_rate=0.5
max_loop=10
li_x=fun_grad(init,learn_rate,max_loop)
# 画图
x = np.arange(0,2*np.pi,0.5)
y=np.cos(x)
plt.scatter(li_x,fun(li_x),c="r")
plt.plot(x,y)
print("init: %d"%init)
print('arg min f(x) of x =', li_x[-1])
print('f(x) =', fun(li_x[-1]))

4. 基于梯度的优化方法–牛顿迭代法
   作业：使用牛顿迭代法求函数 𝑦=𝑐𝑜𝑠(𝑥) 在区间 𝑥∈[0,2𝜋] 的极小值点。
   关于这个方法我有一个问题：牛顿迭代的收敛方向是极值，就算控制了向极小值方向迭代，如果初始点是小于0.5𝜋，大于1.5𝜋，很容易迭代到小于0，大于2𝜋的位置。
   我个人理解是结合二分法进行计算，先找到与极值相近的位置，再迭代到真正的极值点上，如果有更方便的限制方法，请大佬不吝赐教~

def fun_diff3(x):
    return np.tan(x)
def fun_grad(x,max_loop):
    tem=x
    li_x=[x]
    for i in range(max_loop):
        tem=x-fun_diff3(x)
        if np.cos(x)<np.cos(tem):
            tem=x+fun_diff3(x)
        x=tem
        li_x.append(x)
    return li_x
init=6
max_loop=10
li_x=fun_grad(init,max_loop)
# 画图
x = np.arange(0,10*np.pi,0.5)
y=np.cos(x)
plt.scatter(li_x,fun(li_x),c="r")
plt.plot(x,y)
print("init: %d"%init)
print('arg min f(x) of x =', li_x[-1])
print('f(x) =', fun(li_x[-1]))

线性代数复习

1. 从二维推广多维的线性无关及线性相关的含义
   
   具体如上。极大线性无关组是一组向量之间都不相关的方程组，运用极大线性无关组能够有效提取特征信息。
2. 行列式，特征值，特征向量，迹的内在联系
   $$Ax=b$$
   相当于向量x经过A变换得到向量b，
   而A的行列式是指变换的程度，特征值是衡量在特征方向的变换幅度，特征向量则表示特征方向。负数表示变换方向相反。
   迹则是变换的速度，具体对迹的理解可见知乎矩阵的迹的意义
3. 混淆点
   相似变换：
   T为可逆矩阵的情况下
   $$T^{-1}AT=b$$
   二次型：T为向量
   $$q(x)=T^{-1}AT$$

概率论

两个定律比较重要

1. 大数定律
   样本之间独立同分布，样本量大时样本均值收敛于总体均值
2. 中心极限定律
   样本量足够大时，样本均值的分布会趋于正态分布
3. 利用概率估计Π的大小
   随机过程的比较有意思，放两个作业：

import numpy as np
# n 点数
def pi_estimate(n):
    n_rand_X =np.random.uniform(-1,1,n)
    n_rand_Y =np.random.uniform(-1,1,n)
    dis=np.sqrt(n_rand_X**2+n_rand_Y**2)
    n_in=float(len(dis[dis<=1]))
    return 4*n_in/n
a=pi_estimate(1000000)
a
3.140656

4. 利用抽样估计联合分布函数下的概率

# 电子元件寿命服从指数分布，如果测试c个产品，t个寿命超过h的概率
# n->做n次实验
def ele_life(n,c,h,t,lamb):
    times = 0.0
    for i in range(n):
    	# exponential 指数型函数抽样
        c_rand = np.random.exponential(1/lamb,c)
        c_rand_t = len(c_rand[c_rand>h])
        if c_rand_t > t:
            times = times + 1
    return times / n  
print(ele_life(10000,1000,0.6,900,0.1))
print(ele_life(10000,1000,0.6,900,0.2))

数量统计（极大似然估计）

可见西瓜书的推导~
吃瓜教程
除了数学方法推导公式根据公式优化结果
梯度下降的方法，牛顿法都可以使用

随机过程和MCMC

基础的概念

1. 随机变量：随机变量$X$是定义在样本空间$\Omega$上的函数，当x是$X$的观测值时，存在$\Omega$中的$w$使得$x=X(w)$

2. 随机向量：随机向量$(X_1,X_2,...,X_n)$是定义在样本空间$\Omega$上的n元函数，当$(x_1,x_2,...,x_n)$是$(X_1,X_2,...,X_n)$的观测值时，存在w使得$(x_1,x_2,...,x_n)=(X_1(w),X_2(w),...,X_n(w))$，这时称$(x_1,x_2,...,x_n)$为$(X_1,X_2,...,X_n)$的一次观测或者一次实现。

3. 随机过程：设T为$(-\infty ,+\infty )$的子集，若对每个$t\in T$，$X_t$是随机变量，则称随机变量的集合 { X t ∣ t ∈ T } \{X_t|t\in T \} {Xt​∣t∈T}是随机过程。当每个t都有一次观测，那么会形成一条曲线，则称这条曲线为一条轨道或一条轨迹。

4. 有限维分布：对于任何正整数m和T中互不相同的$t_1,...,t_m$，称$(X_{t_1},...,X_{t_m})$的联合分布为随机过程 { X t ∣ t ∈ T } \{X_t|t\in T \} {Xt​∣t∈T}的一个有限维分布，称全体的有限维分布为该随机过程的概率分布。

5. 随机过程的同分布：如果随机过程 { X t ∣ t ∈ T } \{X_t|t\in T \} {Xt​∣t∈T}与随机过程 { Y t ∣ t ∈ T } \{Y_t|t\in T \} {Yt​∣t∈T}有相同的有限维分布，则称他们同分布。

6. 随机过程的独立：如果随机过程 { X t ∣ t ∈ T } \{X_t|t\in T \} {Xt​∣t∈T}中任意选出来的$(X_{t_1},...,X_{t_i})$与从 { Y t ∣ t ∈ T } \{Y_t|t\in T \} {Yt​∣t∈T}中选出来的$(Y_{s_1},...,Y_{s_j})$是相互独立的，则称两个随机过程独立。

7. 随机序列：如果时间集合T是整数，就是一个随机序列（时间序列），记作$X_n$

马尔可夫

马尔可夫性：现在的状态和过去无关，之和上一次有关
一个例子：马尔可夫链

作业

给定下述Rosenbrock函数
$$f(x)=(a-x_1{)}^2+b(x_2-x_1^2{)}^2$$其中
$$x=(x_1,x_2{)}^T\in R^2,a,b\in R$$。试编写程序完成下述工作：
1） 为不同的a,b取值，绘制该函数的3D表面。请问 a,b取值对该表面形状有大的影响吗？，所谓大影响就是形状不再相似。对a,b的取值区间，能否大致给出一个分类，

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import time
%matplotlib inline
from mpl_toolkits.mplot3d import Axes3D
x1 = np.arange(-100, 100, 0.1)
x2 = np.arange(-100, 100, 0.1)
x1, x2 = np.meshgrid(x1, x2)
a = 1
b = 1
# b = 0
# b = -1
z = np.square(a - x1) + b * np.square(x2 - np.square(x1))
fig4 = plt.figure()
ax4 = plt.axes(projection='3d')
ax4.plot_surface(x1, x2, z, alpha=0.3, cmap='winter')  # 生成表面， alpha 用于控制透明度
# 设定显示Lsbel
ax4.set_xlabel('X')
ax4.set_ylabel('Y')
ax4.set_zlabel('Z')
plt.show()

\	$b\in (-\infty ,0)$	$b=0$	$b\in (0,\infty )$
$a\in R$

2） 编写一个算法来找到它的全局最小值及相应的最小解，并在3D图中标出。分析一下你的算法时空效率、给出运行时间。
作业增加了使用梯度下降随机优化的方法，发现从时间效率来说梯度更优秀，但是，从结果上来，牛顿法更准确~

class Rosenbrock():
    def __init__(self):
        self.x1 = np.arange(-100, 100, 0.01)
        self.x2 = np.arange(-100, 100, 0.01)
        #self.x1, self.x2 = np.meshgrid(self.x1, self.x2)
        self.a = 1
        self.b = 1
        self.algorithm_times = 1000
        self.step=1000
        self.answers = []
        self.min_answer_z = []
    # 准备数据
    def data(self):
        z = np.square(self.a - self.x1) + self.b * np.square(self.x2 - np.square(self.x1))
        #print(z.shape)
        return z

    # 随机牛顿
    def snt(self,x1,x2,z,alpha):
        rand_init = np.random.randint(0,z.shape[0])
        x1_init,x2_init,z_init = x1[rand_init],x2[rand_init],z[rand_init]
        x_0 =np.array([x1_init,x2_init]).reshape((-1,1))
        #print(x_0)
        for i in range(self.algorithm_times):
            x_i = x_0 - np.matmul(np.linalg.inv(np.array([[12*x2_init**2-4*x2_init+2,-4*x1_init],[-4*x1_init,2]])),np.array([4*x1_init**3-4*x1_init*x2_init+2*x1_init-2,-2*x1_init**2+2*x2_init]).reshape((-1,1)))
            x_0 = x_i
            x1_init = x_0[0,0]
            x2_init = x_0[1,0]
        answer = x_0
        return answer
    # 梯度下降
    def grad(self,x1,x2,z,alpha):
        rand_init = np.random.randint(0,z.shape[0])
        x1_init,x2_init,z_init = x1[rand_init],x2[rand_init],z[rand_init]
        x_0 =np.array([x1_init,x2_init]).reshape((-1,1))
        #print(x_0)
        for i in range(self.algorithm_times):            
            x_i = x_0 - alpha*np.array([4*x1_init**3-4*x1_init*x2_init+2*x1_init-2,-2*x1_init**2+2*x2_init]).reshape((-1,1))
            dis=x_0-x_i
            if np.all(dis<=np.array([1e-4,1e-4]).reshape((-1,1))):
                break
            x_0 = x_i
            x1_init = x_0[0,0]
            x2_init = x_0[1,0]            
        answer = x_0
        return answer
    # 绘图
    def plot_data(self,min_x1,min_x2,min_z):
        x1 = np.arange(-100, 100, 1)
        x2 = np.arange(-100, 100, 1)
        x1, x2 = np.meshgrid(x1, x2)
        a = self.a
        b = self.b
        z = np.square(a - x1) + b * np.square(x2 - np.square(x1))
        fig4 = plt.figure()
        ax4 = plt.axes(projection='3d')
        ax4.plot_surface(x1, x2, z, alpha=0.3, cmap='winter')  # 生成表面， alpha 用于控制透明度
        ax4.scatter(min_x1,min_x2,min_z,c='r')
        ax4.text(min_x1+10,min_x2+10,min_z,r'({:.2f},{:.2f},{:.2f})'.format(min_x1,min_x2,min_z),fontsize=10)
        # 设定显示范围
        ax4.set_xlabel('X')
        ax4.set_ylabel('Y')
        ax4.set_zlabel('Z')
        plt.show()
    # 开始
    def start(self):
        times = int(input("请输入需要随机优化的次数："))
        alpha = float(input("请输入随机优化的步长"))
        z = self.data()
        start_time = time.time()
        for i in range(times):
            answer = self.snt(self.x1,self.x2,z,alpha)
#             answer = self.grad(self.x1,self.x2,z,alpha)
            self.answers.append(answer)
        min_answer = np.array(self.answers)
        for i in range(times):
            self.min_answer_z.append((1-min_answer[i,0,0])**2+(min_answer[i,1,0]-min_answer[i,0,0]**2)**2)
        optimal_z = np.min(np.array(self.min_answer_z))
        optimal_z_index = np.argmin(np.array(self.min_answer_z))
        optimal_x1,optimal_x2 = min_answer[optimal_z_index,0,0],min_answer[optimal_z_index,1,0]
        end_time = time.time()
        running_time = end_time-start_time
        print("优化的时间:%.2f秒!" % running_time)
        self.plot_data(optimal_x1,optimal_x2,optimal_z)
if __name__ == '__main__':
    snt = Rosenbrock()
    snt.start()

牛顿法结果	梯度下降结果