集成学习学习笔记——数学基础（1）

一、高等数学与线性代数

1.Jacobi矩阵

假设F：Rn——>Rm，则F的Jacobi矩阵为【yi/xj】(m*n)

2.Hessian矩阵         

假设F：Rn——>R，则F的Hessian矩阵为

二、概率论与随机过程初步

概率空间表示为,其中Omega为全体trial结果的集合，F为Omega的幂集的子集（对事件集合的），P为概率测度。随机变量则为Borel可测函数

随机过程略

三、拒绝采样与MCMC采样

1.拒绝采样

当有一个复杂概率分布，我们需要依照该分布进行随机抽样，以下方法称为拒绝采样：

Step1：取易于抽样的分布（例如均匀分布、正态分布等），且

Step2：While 取样数量未达到目标：

             1）从中取样得到

             2）依照分布抽样得到

             3）若，则说明落在了f的分布范围内，接受此样本，否则拒绝此样本

2.Markov链

1）Markov过程

2）Markov链收敛性质

对非周期Markov链有转移矩阵P，且任意状态间连通（不可约，整个Markov图中只有一个连通类），有以下性质，其中表示状态空间的分布：

【I】

【II】是方程的唯一非负解

3）基于Markov链采样

如果得到了某个平稳分布所对应的Markov链转移矩阵，我们就很容易采集到这个平稳分布的样本集，假设初始状态的概率分布为，第一轮转移后为，。。。，第i轮后为，直至第n轮后达到平稳，即：

                                

采样过程：

Step1：输入Markov状态转移矩阵P，设定转移次数阈值，采样样本数

Step2：从简单分布（uniform,normal）采样得到初始状态值

Step3：loop t=0 to +-1：从条件概率分布中采样得到

样本集合即为平稳分布对应的样本集

3.MCMC

【Markov链细致平稳条件】对非周期Markov链的状态转移矩阵P和概率分布，，则称概率分布是状态转移矩阵P的平稳分布。

【MCMC采样】

一般情况下，目标平稳分布和某一个Markov状态转移矩阵Q不满足细致平稳条件，即：

因此做如下改造：

称为接受率，处于[0,1]，可以理解为一个概率值，目标矩阵P可以通过任意一个Markov状态转移矩阵Q以特定的接受率获得。

Algo：

1）输入任意选定的Markov链状态转移矩阵Q，平稳分布，设定转移次数阈值，采样样本数

2）从简单分布（uniform,normal）采样得到初始状态值

3）loop t=0 to +-1：
        a）从条件概率分布中采样得到样本

        b）从均匀分布中采样

        c）如果，则接收转移，即

        d）否则

样本集合即为平稳分布对应的样本集

Disadvantages：

3）c）中，接受率一般十分小，导致大部分转移都被拒绝，马尔可夫链不收敛，需要取得很大

【Metropolis-Hastings采样】

对接受率做等比缩放即可，例如

 问题一般是太小了，但以上等式两侧同时扩大n倍依旧成立，故对接受率做如下调整：

Algo：

1）输入任意选定的Markov链状态转移矩阵Q，平稳分布，设定转移次数阈值，采样样本数

2）从简单分布（uniform,normal）采样得到初始状态值

3）loop t=0 to +-1：
        a）从条件概率分布中采样得到样本

        b）从均匀分布中采样

        c）如果，则接收转移，即

        d）否则

样本集合即为平稳分布对应的样本集

Notation：选择对称的Markov状态转移矩阵Q，接受率可转换为

四、Exercise

给定下述Rosenbrock函数，，其中，。试编写程序完成下述工作：

1. 为不同的a,b取值，绘制该函数的3D表面。请问 a,b取值对该表面形状有大的影响吗？，所谓大影响就是形状不再相似。对a,b的取值区间，能否大致给出一个分类，像下面这样给出一张表：

b>0	b<0	b=0

      2.编写一个算法来找到它的全局最小值及相应的最小解，并在3D图中标出。分析一下你的算法            时空效率、给出运行时间。

思路：利用梯度下降法进行搜索

初始化参数a=1，b=1，代码如下

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

plt.style.use("ggplot")


def rosenbrock_func(param_1, param_2, a, b):
    return (a - param_1) ** 2 + b * (param_2 - param_1 ** 2) ** 2


def rosenbrock_delta(param_1, param_2, a, b, param_index):
    if param_index == 1:
        return 2 * (param_1 - a) + 4 * b * (param_1 ** 2 - param_2) * param_1
    else:
        return 2 * b * (param_2 - param_1 ** 2)


def gradient_descent(x1_primal, x2_primal, a, b, learning_rate=0.1, error_threshold=1e-5, max_iter=1e5):
    iter_num, error, y_old = 0, float('inf'), rosenbrock_func(x1_primal, x2_primal, a, b)
    gd_x1, gd_x2, gd_y, error_list = [x1_primal], [x2_primal], [y_old], [error]
    while error > error_threshold and iter_num < max_iter:
        iter_num += 1
        x1_primal = x1_primal - learning_rate * rosenbrock_delta(x1_primal, x2_primal, a, b, 1)
        x2_primal = x2_primal - learning_rate * rosenbrock_delta(x1_primal, x2_primal, a, b, 2)
        y_new = rosenbrock_func(x1_primal, x2_primal, a, b)
        error = abs(y_new - y_old)
        y_old = y_new
        gd_x1.append(x1_primal)
        gd_x2.append(x2_primal)
        gd_y.append(y_new)
        error_list.append(error)
    return [gd_x1, gd_x2, gd_y, iter_num, error_list]


def func_plot(x1, x2, y, ex_x1=None, ex_x2=None, ex_y=None):
    fig = plt.figure()
    ax1 = plt.axes(projection='3d')
    ax1.plot_surface(x1, x2, y, alpha=0.3, cmap='rainbow')
    if ex_x1 and ex_x2 and ex_y:
        ax1.scatter3D(ex_x1, ex_x2, ex_y, c='red', norm=1, )
    ax1.set_xlabel('X1')
    ax1.set_ylabel('X2')
    ax1.set_zlabel('Y')
    plt.show()
    return


if __name__ == '__main__':
    param1_array, param2_array = np.arange(-10, 10, 0.01), np.arange(-10, 10, 0.01)
    x1_list, x2_list = np.meshgrid(param1_array, param2_array)
    y_list = rosenbrock_func(x1_list, x2_list, 1, 1)
    gd_x1, gd_x2, gd_y, iter_num, error_list = gradient_descent(2, 2, 1, 1)
    print(gd_x1[-1], gd_x2[-1], gd_y[-1], iter_num, error_list[-1])
    func_plot(x1_list, x2_list, y_list, gd_x1[-1], gd_x2[-1], gd_y[-1])

迭代次数53次，求解得到全局最小值0.0001108，x1=1.00978，x2=1.02354

在函数中表现为：