精华内容
下载资源
问答
  • Actor-Critic核心在Actor以下分三个部分介绍Actor-Critic方法,分别为(1)基本的Actor算法(2)减小Actor的方差 (3)Actor-Critic。仅需要强化学习的基本理论和一点点数学知识。基本的Actor算法Actor基于策略梯度,策略...

    Actor-Critic核心在Actor

    以下分三个部分介绍Actor-Critic方法,分别为(1)基本的Actor算法(2)减小Actor的方差 (3)Actor-Critic。仅需要强化学习的基本理论和一点点数学知识。

    基本的Actor算法

    Actor基于策略梯度,策略被参数化为神经网络,用

    表示。

    迭代的方向是最大化周期奖励的期望,目标函数表示为:

    其中

    代表一个采样周期,

    代表序列出现的概率。求

    的梯度可得:

    其实仅是将

    写成了

    而已,目的是期望仍然为对采样序列的期望。

    表示序列的奖励和。

    那么关键的地方来了。

    代表整个序列出现的概率,整个序列的出现包括两个影响因素:策略

    和由环境决定的状态转移

    写成下式:

    很好理解,”初始状态的概率“乘以很多次(”动作选择概率“乘以”由环境决定的状态转移概率“),后两项需要连乘,直到周期结束。根据上面的式子,就可以求

    中出现的梯度项

    的值。在求

    时,上式中的

    由环境决定,与我们关注的策略变量

    无关(就是不管你采用什么策略,环境照样那么转移),因此在计算中被略去(如果这几项无法略去,那么在无环境(model-free)的情形下,将无法计算关于策略的梯度)。推导可得策略的梯度,如下:

    虽然式子看起来比较复杂,其实很好理解。(1)期望仍是对整个序列出现的期望,在强化学习中通过采样进行估计即可,例如Monte Carlo采样。(2)期望中第一项是

    取对数再求导后的结果,有关环境因素

    被消去,仅剩与策略

    相关的因素。(3)期望第二项是

    的展开,此处为了书写方便,仅使用所有奖励的和,在实际问题中可以使用average reward、discount reward等。

    以上公式与Monte Carlo算法结合起来便构成了”REINFORCE策略梯度“算法。该算法分为三个步骤:

    (1)采样

    个序列

    (只有采样才能估计梯度的期望,采样多了就估计准了)

    (2)计算策略梯度:

    。就是对每一个采样序列分别计算梯度随后取平均。

    (3)更新策略的参数:

    是学习率,由于是最大化

    ,因此是梯度上升。

    重复(1)(2)(3)三个步骤,直到收敛。可以看出,该算法是on-policy的。

    以上,便构造了一个基本的策略梯度算法。

    在实现中,需要根据实际情况设计策略网络来表示

    ,在给定状态下,输出动作选择的概率。分别讨论离散动作空间和连续动作空间下的网络设计。

    (1)在离散动作空间中,输入为状态的表示,输出节点与动作个数相等,后接Softmax层。

    (2)在连续动作空间中,输入为状态的表示,输出的设计方式有多种。一般假设每个动作的输出服从高斯分布,因此可以输出每个动作的均值。动作之间可以共用方差或各自分别学习方差。近期也有研究指出输出使用Beta分布比Gussian分布要好。

    减小Actor的方差

    Causality方法

    策略梯度的公式变为

    即后一项关于的奖励的累加只累加当前时间步

    之后的奖励,直观的理解是

    时刻策略的改变不会影响到

    之前时间步的奖励。在上式中,对

    使用的是Monte Carlo估计方法,这种方法方差大。将该项可以写成

    ,Q值代表未来的累计奖励的期望,可以使用值函数近似的方法来估计

    ,从而进一步减少方差。这是Actor-Critic的一种。此时策略梯度的公式变为:

    Baseline方法

    如果希望在上式的基础上,进一步减少方差,那么可以为

    添加baseline,将baseline记为

    ,则策略梯度的公式变为:

    可以证明,只有在

    与动作

    无关的情况下,上述改进才与之前的策略梯度公式等价。

    一般选择为状态

    的值函数,即

    Off-policy

    REINFORCE算法是一个on-policy算法,每次改变策略时都需要重新采集样本,因此样本利用效率低。使用importance Sampling方法可以将其转变为off-policy,这里不再赘述,转变后策略梯度的公式变为:

    其中,

    是采样策略的参数。importance sampling项的连乘范围从周期的开始一直到当前时间步

    Actor-Critic

    我们在baseline的基础上继续进行分析。当baseline是状态值函数时,策略梯度可以写成

    其中,

    称为Advantage Function。而

    如果不考虑状态转移概率,用采样的方式来估计状态转移,则在当前策略参数下,

    因此策略梯度公式可以进一步写成

    在上式中,我们需要估计状态值函数

    的值。用于估计

    的部分被称为Critic.

    熟悉Q学习的童鞋应该知道,估计状态值函数的方法有很多种,在此不再赘述。

    Actor-Critic的基本流程为:

    采样

    更新Critic参数

    根据Critic计算Advantage Function

    更新Actor参数

    Actor-Critic方法也存在问题。虽然减少了Actor在计算策略梯度中的方差,但是由于Critic开始时对V函数估计不准,会导致Actor对策略梯度的估计存在偏差。特别是在On-line更新中,仅通过单个样本来估计状态值函数是不现实的,因此有人提出了改进,用多个并行的Actor-Critic进行学习,就是大名鼎鼎的A3C算法。另外还有人对baseline进行研究,提出动作状态值函数也可以直接作为baseline,不过需要添加补偿项(Q-prop)。DDPG也是一种Actor-Critic算法,其特点为策略参数的更新的目标不是

    ,而是最大化

    值,使Q值最大化的策略可以使用确定性的策略,因此成为确定性策略梯度。

    在网络设计方面,Actor和Critic可以使用不同的网络,这样学习稳定。也可以使用同一个网络,共享底层的特征。

    以上。

    如有疏漏之处,欢迎指正。

    展开全文
  • 来自|PaperWeekly作者|邓妍蕾学校|...日常中我们会采样的方法采集样本,进行近似的数值计算,比如计算样本均值,方差期望。虽然许多常见的概率密度函数(t 分布,均匀分布,正态分布),我们都可以直接在 ...

    来自|PaperWeekly

    作者|邓妍蕾

    学校|香港大学硕士

    研究方向|NLP、语音识别

    概览

    马尔科夫蒙特卡洛法(Markov Chain Monte Carlo, MCMC)经常用在贝叶斯概率模型的推理和学习中,主要是为了解决计算困难的问题。

    日常中我们会用采样的方法采集样本,进行近似的数值计算,比如计算样本均值,方差,期望。虽然许多常见的概率密度函数(t 分布,均匀分布,正态分布),我们都可以直接在 Numpy,Scikit-Learn 中找到,但很多时候,一些概率密度函数,特别是高维的概率密度函数,并不是常见的分布,这个时候我们就需要用到 MCMC。

    在开始马尔科夫蒙特卡洛法之前,我们先简单介绍一些蒙特卡洛法和马尔科夫链。

    蒙特卡洛法 Monte Carlo

    蒙特卡洛法是比较宽泛的一系列算法的统称(想了解可自行 google),它的特点是假设概率分布已知,通过重复的随机采样来获得数值结果。比如根据大数定理,我们可以用采样得到的样本计算得到的样本均值来估计总体期望(例子 1)。又比如,积分的运算往往可以表示为随机变量在一个概率密度函数分布上的期望(例子 2)。

    例子 1:假设有随机变量 x,定义域 X,其概率密度函数为 p(x),f(x) 为定义在 X 上的函数,目标是求函数 f(x) 关于密度函数 p(x) 的数学期望

    。蒙特卡洛法根据概率分布 p(x) 独立地抽样 n 个样本
    ,得到近似的 f(x) 期望为:

    40139d5ec6ec1f57c2fa1574cedf7261.png

    例子 2: 假设我们想要求解 h(x) 在 X 上的积分:

    5bf0341b27fd736c16305ba0fd7c227c.png

    我们将 h(x) 分解成一个函数 f(x) 和一个概率密度函数 p(x) 的乘积,进而又将问题转换为求解函数 f(x) 关于密度函数 p(x) 的数学期望

    b0fc7067f57d9736d660564895bc7736.png

    这里,f(x) 表示为

    ,则有:

    1a27ed02b9e24cbca27b8a4d083c36e4.png

    更一般的,假设我们想要求解

    ,熟悉积分的同学肯定已经知道答案为
    ,那么如何用采样的方法来得到这个值呢?

    ,0<x<10,那么
    import random
    numSamples = 10000
    # 在 0-10的均匀分布内采样
    samples = [random.uniform(0,10) for _ in range(num_samples)]
    f_samples = [10 * sample ** 2 for sample in samples]
    result = 1/10000.0 * sum(f_samples)
    #>>> result
    #333.7766822849899

    对于复杂的 h(x),这种方法计算起来显然就更加方便了(特别是忘记积分怎么算的同学)。

    到这里为止,我们简单的介绍了蒙特卡洛方法,但是依旧没有提到要怎么利用复杂的概率密度函数进行采样。接下来我们来看一下接受-拒绝法(accept-reject sampling method),它也是蒙特卡洛法中的一种类型适用于不能直接抽样的情况。

    接受-拒绝法

    假设有一个非常复杂不常见的分布 p(x),我们没有现成的工具包可以调用后进行采样,那么我们可以用我们已经有的采样分布比如高斯分布作为建议分布(proposal distribution),用 q(x) 表示,来进行采样,再按照一定的方法来拒绝当前的样本,使得最后得到的样本尽可能的逼近于 p(x)。

    首先我们需要找到一个常数 k 使得 kq(x) 一定大于等于 p(x), 也就是如图所示(图摘自 [2]),p(x) 在 q(x) 下面。接着对 q(x) 进行采样,假设得到的样本为

    。然后我们按照均匀分布在
    中采样得到 u。如果 u 落到了图中的灰色区域,则拒绝这次采样,否则接受样本
    。重复这个过程得到一系列的样本

    5039f07e515d8d2ab4ae4ef45a411e7a.png

    在这个过程中,我们可以发现只有当建议分布 q(x) 和复杂分布 p(x) 重合的尽可能多的地方的样本更有可能被接受。那么相反的,如果他们重合部分非常少,那么就会导致拒绝的比例很高,抽样的效率就降低了。很多时候我们时候在高维空间进行采样,所以即使 p(x) 和 q(x) 很接近,两者的差异也可能很大。

    我们可以发现接受-拒绝法需要我们提出一个建议分布和常量,而且采样的效率也不高,那么我们就需要一个更一般的方法,就是我们的马尔科夫蒙特卡洛法。不过 MCMC 依旧有个问题, 它的抽样样本不是独立的。到了 Gibbs Sampling 的部分,我们可以看到它做了一些小 trick,来假装使得样本独立。

    马尔科夫链马尔科夫链的定义和概念相信大家都很熟悉了(不熟悉的请 google)。在这里,我们回顾一下几个重要的知识点:

    平稳分布

    设有马尔科夫链

    , 其状态空间为 S,转移概率矩阵为
    ,如果存在状态空间 S 上的一个分布:

    97dd3dc5cd8fb47a1ea297e596486acd.png

    使得 π=Pπ,则称 π 为马尔科夫链

    的平稳分布。

    我们可以理解为,当一个初始分布为该马尔科夫链的平稳分布的时候,接下来的任何转移操作之后的结果分布依然是这个平稳分布。注意,马尔科夫链可能存在唯一平稳分布,无穷多个平稳分布,或者不存在平稳分布。

    其它定理:

    1. 不可约且非周期的有限状态马尔科夫链,有唯一平稳分布存在。

    2. 不可约、非周期且正常返的马尔科夫链,有唯一平稳分布存在。

    其中,不可约和正常返大家请自行查阅相关定义。直观上可以理解为,任何两个状态都是连通的,即从任意一个状态跳转到其它任意状态的概率值大于零。

    遍历定理

    设有马尔科夫链

    , 其状态空间为 S,若马尔科夫链 X 是不可约、非周期且正常返的,则该马尔科夫链有唯一平稳分布
    ,且转移概率的极限分布是马尔科夫链的平稳分布。

    3d5ead4b2852ae5684e6536278319295.png

    直观上,P 最后会收敛成这个样子:

    c1782e323fbe57e5a2aaf564acb2afdb.png

    注意:这里的

    表示的是第 j 个状态转移到第 i 个状态的概率。也就是说每一列表示这个状态转移到其它状态的概率。我们可以理解为,当时间趋于无穷时,马尔科夫链的状态分布趋近于平稳分布。另外,无论初始的分布 π′ 是什么,经过了 n 次转移后,即
    ,最后都会收敛到它的平稳分布 π(原理待补充)。结合这两点,我们就可以用马尔科夫链进行采样。

    首先我们随机一个样本

    ,基于条件概率(转移概率)
    采样
    ,因为我们需要转移一定次数后才会收敛到我们的平稳分布,所以比如我们设定了 m 次迭代后为平稳分布,那么
    即为平稳分布对应的样本集。

    但是,要怎么确定平稳分布 π(我们希望采样的复杂分布)的马尔科夫链状态转移矩阵或者转移核 P 呢?在开始 MCMC 采样之前我们还需要回顾两个知识点: 可逆马尔科夫链和平衡细致方程。

    可逆马尔科夫链

    设有马尔科夫链

    ,其状态空间为 S,转移概率矩阵为 P,如果有状态分布
    ,对于任意状态 i,j∈S,对于任意一个时刻 t 满足:

    552a468c12a094f86eaf318d2cedcbf3.png

    或简写为:

    70614030564bb8eae7361999f9150e37.png

    该式也被称之为细致平衡方程。 定理:满足细致平衡方程的状态分布 π 就是该马尔科夫链的平稳分布,即 Pπ=π。因此,可逆马尔科夫链一定有唯一平稳分布,所以可逆马尔科夫链满足遍历定理的条件。

    马尔科夫链蒙特卡洛法

    我们先来看一下 MCMC 方法的大致思路:在随机变量 x 的状态空间 S 上定义一个满足遍历定理的马尔科夫链

    ,使其平稳分布就是抽样的目标分布 p(x)。然后在这个马尔科夫链上进行随机游走,每个时刻得到一个样本。根据遍历定理,当时间趋向于无穷时,样本的分布趋近于平稳分布,样本的函数均值趋近函数的数学期望。

    读者可能会和笔者刚读到这段话的时候也有些疑惑:1)如何定义这个能满足条件的马尔科夫链呢?2)这里的随机游走(转移核)采样具体是怎么实现的?

    对于第一个问题,定义满足条件的马尔科夫链,大家也许猜到了可以用可逆马尔科夫链的性质定义一些特殊的转移核来保证遍历定理成立。但是我们又该如何保证这个特殊的转移核最后可以走到我们期望采样的复杂分布 P(x) 呢?看起来还是一头雾水。我们接着往下看。

    Metropolis-Hastings

    Metropilis-Hastings 算法是马尔科夫链蒙特卡洛法的代表算法。假设要抽样的终极概率分布是 p(x),它采用的转移核为:

    63fdd47b9bd133010aed8cfd1cddafae.png

    其中 q(x,x’) 是另一个卡尔科夫链的转移核,并且是一个容易抽样的分布,被称之为建议分布。而 α(x,x′) 被称为接受分布:

    f5364c39d87743855f5f89212cdffd13.png

    到这里为止,读者可能觉得好像又和之前的拒绝-接受法很相似了呢。但是不同的地方是之前我们在拒绝-接受法里需要定义的建议分布还是比较难设定的,因为它需要满足一定的条件才可以。但是这里的建议分布相是一个比较容易抽样的分布。同理,我们根据建议分布 q(x,x′) 来进行随机游走,产生样本后,按照接受分布 α(x,x′) 来确定是否要进行转移。

    那接下来,我们需要解决的主要问题是,如何证明通过这个方式最后生成的样本是符合 p(x) 分布的呢?也就是说,如何证明这个转移核 p(x,x′) 满足遍历定理以及最后的平稳分布是 p(x) 呢?

    证明略。其实我们只要能证明这个构造出来的转移核 p(x,x′) 对应的马尔科夫链是可逆的,并且其对应的平稳分布就是 π 即可。也就是说需要证明:

    f52928c4ff2b2f416081d54c0d044987.png

    证明:

    8fc50903543427f1e04c8346aa796b82.png

    并且根据细致平衡方程,p(x) 即为这个转移核 p(x,x′) 的平衡分布。

    在给出 M-H 方法的总体流程之前,我们先来看几个特殊的建议分布。

    建议分布

    对称建议分布:假设我们的建议分布是对称的,即 q(x,x′)=q(x′,x),那么我们的接受分布 α(x,x′) 可以写成:

    a14ba6d582cb27262782b67b55e56653.png

    特别地,q(x,x′)=q(|x−x′|) 被称为随机游走 Metropolis 算法,例子:

    121dfbca3b4600d08db96c64c73873d3.png

    读者也很容易发现,当正态分布的方差为常数,均值为 x,参数为 x′ 的这些转移核都满足这种类型。这种类型的转移核的特点是,当 x′ 在均值 x 附近的时候,其概率也就越高。

    独立抽样:前面的抽样过程中,可以发现下一个样本总是依赖于前一个样本。那么我们假设设定的建议分布 q(x,x’) 与当前状态 x 无关,即 q(x,x′)=q(x′), 此时的接受分布 α(x,x′) 可以写成:

    a063298710778a08effb16a3d519c78d.png

    书上说,这样的抽样虽然简单,但是收敛速度慢,通常选择接近目标分布 p(x) 的分布作为建议分布 q(x)。

    接下来,我们看一下 M-H 方法的总体过程。

    Metropolis-Hastings 算法步骤

    输入:任意选定的建议分布(状态转移核)q,抽样的目标分布密度函数 p(x),收敛步数 m,需要样本数 n。

    Step 1:随机选择一个初始值

    ,样本集合samples=[]

    Step 2:for t=0 to m+n:

    • 按照建议分布 q(x,x′) 随机抽取一个候选状态 x′
    • 计算接受概率:

    d895e4e80044484b392227c0f9625c6a.png
    • 从区间 (0,1) 中按均匀分布随机抽取一个数 u。若 u≤α(x,x′),则状态 x=x′,否则,x 保持不变
    • if (t >= m) 将 x 加入到 samples 中

    Step 3:返回样本集合 samples。一个简单的例子:

    # -*- coding:utf-8 -*-
    
    import random
    import numpy as np
    import matplotlib.pyplot as plt
    
    
    def mh(q, p, m, n):
        # randomize a number
        x = random.uniform(0.1, 1)
        for t in range(0, m+n):
            x_sample = q.sample(x)
            try:
                accept_prob = min(1, p.prob(x_sample)*q.prob(x_sample, x)/(p.prob(x)*q.prob(x, x_sample)))
            except:
                accept_prob = 0
    
            u = random.uniform(0, 1)
    
            if u < accept_prob:
                x = x_sample
    
            if t >= m:
                yield x
    
    
    class Exponential(object):
        def __init__(self, scale):
            self.scale = scale
            self.lam = 1.0 / scale
    
        def prob(self, x):
            if x <= 0:
                raise Exception("The sample shouldn't be less than zero")
    
            result = self.lam * np.exp(-x * self.lam)
            return result
    
        def sample(self, num):
            sample = np.random.exponential(self.scale, num)
            return sample
    
    
    # 假设我们的目标概率密度函数p1(x)是指数概率密度函数
    scale = 5
    p1 = Exponential(scale)
    
    
    class Norm():
        def __init__(self, mean, std):
            self.mean = mean
            self.sigma = std
    
        def prob(self, x):
            return np.exp(-(x - self.mean) ** 2 / (2 * self.sigma ** 2.0)) * 1.0 / (np.sqrt(2 * np.pi) * self.sigma)
    
        def sample(self, num):
            sample = np.random.normal(self.mean, self.sigma, size=num)
            return sample
    
    # 假设我们的目标概率密度函数p1(x)是均值方差分别为3,2的正态分布
    p2 = Norm(3, 2)
    
    
    class Transition():
        def __init__(self, sigma):
            self.sigma = sigma
    
        def sample(self, cur_mean):
            cur_sample = np.random.normal(cur_mean, scale=self.sigma, size=1)[0]
            return cur_sample
    
        def prob(self, mean, x):
            return np.exp(-(x-mean)**2/(2*self.sigma**2.0)) * 1.0/(np.sqrt(2 * np.pi)*self.sigma)
    
    
    # 假设我们的转移核方差为10的正态分布
    q = Transition(10)
    
    m = 100
    n = 100000 # 采样个数
    
    simulate_samples_p1 = [li for li in mh(q, p1, m, n)]
    
    plt.subplot(2,2,1)
    plt.hist(simulate_samples_p1, 100)
    plt.title("Simulated X ~ Exponential(1/5)")
    
    samples = p1.sample(n)
    plt.subplot(2,2,2)
    plt.hist(samples, 100)
    plt.title("True X ~ Exponential(1/5)")
    
    simulate_samples_p2 = [li for li in mh(q, p2, m, n)]
    plt.subplot(2,2,3)
    plt.hist(simulate_samples_p2, 50)
    plt.title("Simulated X ~ N(3,2)")
    
    
    samples = p2.sample(n)
    plt.subplot(2,2,4)
    plt.hist(samples, 50)
    plt.title("True X ~ N(3,2)")
    
    plt.suptitle("Transition Kernel N(0,10)simulation results")
    plt.show()

    代码运行结果:

    f521a5e6714d0b7c5f9efe9e57f792f7.png

    多元情况

    在很多情况下,我们的目标分布是多元联合概率分布

    ,我们可以用条件分布概率的比来计算联合概率的比,从而提升计算效率。即:

    1935b429958d802770320abcb8176b42.png

    其中:

    e6220bddda4aa28da5ee8b73a87d4787.png

    并且

    被称为满条件分布(full conditional distribution)。

    而在用转移核对多元变量分布进行抽样的时候,可以对多元变量的每一个变量的条件分布一次分别进行抽样,从而实现对整个多元变量的一次抽样。

    针对多元变量的情况,我们来讨论一下其中的一个特例,吉布斯采样(Gibbs Sampling)。

    Gibbs Sampling

    吉布斯采样(Gibbs Sampling)适用于多元变量联合分布的抽样和估计。比较特殊的地方是它的建议分布是当前变量

    的满条件概率分布:

    304f95c3714ff37344135afa5a787c78.png

    这是,接受概率 α=1:

    1403f53800f67a14dc751280d42fbad3.png

    其中,

    因为不难想象,假设对第 j 个变量采样前
    ,对其采完样后
    ,那么在 x′ 中分别去掉第 j 项, 剩下的其余变量全部相同,所以它们的概率值也相同。

    同理:

    ed1d62a654817292d51f6b2817a8893e.png

    我们可以发现因为 α 接受分布永远等于 1,那么吉布斯采样是不会拒绝样本的。

    这里需要注意的是,这里的目标分布的单个元的分布至少应该是可以采样的,当然如果单个元很复杂的话,还是可以再用一次 M-H 的方法来为这个单元的分布进行采样(存疑)。

    吉布斯采样算法

    输入:抽样的目标分布密度函数 p(x),转移核

    收敛步数 m,需要样本数 n。

    Step 1:随机选择一个样本

    ,样本集合 samples=[]

    Step 2:

    for i=1 to m+n:

    c1848e317b347c4b584161969b1113a5.png

    for j=1 to k:

    根据

    抽取
    , 并赋值给

    if (i >= m) 将

    加入到 samples 中。

    Step 3:返回 samples。

    来看个例子吧。假设我们有二维正态分布:

    a023862ce300f09a9f0b22f2f8845b5b.png

    假设我们期望抽样的二元正态分布是 μ=(5,8),协方差矩阵为:

    645d07257ea596b397e75fb24c6d8d45.png
    import matplotlib.pyplot as plt
    from mpl_toolkits.mplot3d import Axes3D
    import numpy as np
    
    
    class Transition():
        def __init__(self, mean, cov):
            self.mean = mean
            self.sigmas = []
            for i in range(K):
                self.sigmas.append(np.sqrt(cov[i][i]))
            self.rho = cov[0][1]/(self.sigmas[0] * self.sigmas[1])
    
        def sample(self, id1, id2_list, x2_list):
            id2 = id2_list[0]  # only consider two dimension
            x2 = x2_list[0]  # only consider two dimension
            cur_mean = self.mean[id1] + self.rho*self.sigmas[id1]/self.sigmas[id2] * (x2-self.mean[id2])
            cur_sigma = (1-self.rho**2) * self.sigmas[id1]**2
            return np.random.normal(cur_mean, scale=cur_sigma, size=1)[0]
    
    
    def gibbs(p, m, n):
        # randomize a number
        x = np.random.rand(K)
        for t in range(0, m+n):
            for j in range(K):
                total_indexes = list(range(K))
                total_indexes.remove(j)
                left_x = x[total_indexes]
                x[j] = p.sample(j, total_indexes, left_x)
    
            if t >= m:
                yield x
    
    
    mean = [5, 8]
    cov = [[1, 0.5], [0.5, 1]]
    K = len(mean)
    q = Transition(mean, cov)
    m = 100
    n = 1000
    
    gib = gibbs(q, m, n)
    
    simulated_samples = []
    
    x_samples = []
    y_samples = []
    for li in gib:
        x_samples.append(li[0])
        y_samples.append(li[1])
    
    
    fig = plt.figure()
    ax = fig.add_subplot(131, projection='3d')
    
    hist, xedges, yedges = np.histogram2d(x_samples, y_samples, bins=100, range=[[0,10],[0,16]])
    xpos, ypos = np.meshgrid(xedges[:-1], yedges[:-1])
    xpos = xpos.ravel()
    ypos = ypos.ravel()
    zpos = 0
    
    dx = xedges[1] - xedges[0]
    dy = yedges[1] - yedges[0]
    dz = hist.flatten()
    
    ax.bar3d(xpos, ypos, zpos, dx, dy, dz, zsort='average')
    
    ax = fig.add_subplot(132)
    ax.hist(x_samples, bins=50)
    ax.set_title("Simulated on dim1")
    
    ax = fig.add_subplot(133)
    ax.hist(y_samples, bins=50)
    ax.set_title("Simulated on dim2")
    plt.show()

    代码运行结果:

    336c99636132fb21de62c01c88e3626d.png

    Tricks

    从马尔科夫链蒙特卡洛法中我们得到的样本序列,相邻点是相关的,因此如果需要独立样本,我们可以在该样本序列中再次进行随机抽样,比如每隔一段时间抽取一次样本,将这样得到的样本集合作为独立样本集合。

    另外,关于“燃烧期”,即过多少步后马尔科夫链是趋向稳定分布的,通常是经验性的。比如书上说,可以游走一段时间后采集的样本均值和过一段时间后采集的样本均值做个比较,如果均值稳定,那可以认为马尔科夫链收敛。

    总结

    本文回顾了马尔科夫链蒙特卡洛法的相关基本知识点、Metropolis-Hastings, Gibbs sampling 基本原理。下一篇文章,我们一起来看看 MCMC 大法在实际应用中的一些例子。

    Reference

    [1] 李航. 统计学习方法第二版[J]. 2019.

    [2] https://www.cnblogs.com/pinard/p/6625739.html

    关于PaperWeekly

    PaperWeekly 是一个推荐、解读、讨论、报道人工智能前沿论文成果的学术平台。如果你研究或从事 AI 领域,欢迎在公众号后台点击「交流群」,小助手将把你带入 PaperWeekly 的交流群里。

    加入社区:http://paperweek.ly

    微信公众号:PaperWeekly

    新浪微博:@PaperWeekly

    展开全文
  • 期望方差线性化的卡尔曼滤波器称作扩展卡 尔曼滤波器(Extended Kalman Filter), 简称EKF。 同泰勒级数类似, 面对非线性关系时, 我们可以通过求过程和量测方 程的偏导来线性化并计算当前估计。 我们将第一...
  • DeepLearning的数学基础

    2018-12-06 11:28:39
    第⼀一章 数学基础 1.1 标量量、向量量、矩阵、张量量之间的联系 1.2 张量量与矩阵的区别? 1.3 矩阵和向量量相乘结果 1.4 向量量和矩阵的范数归纳 1.5 如何判断⼀一个矩阵...1.17 期望、⽅方差、协⽅方差、相关系数总结
  • 统计模拟及其R实现

    2014-11-19 12:39:21
    在对条件期望、条件方差、Poisson过程和Markov链的基本知识进行简单介绍之后,介绍了如何利用计算机产生随机数以及如何利用这些随机数产生任意分布的随机变量、随机过程等知识;介绍了一些分析统计数据的方法和技术...
  • 在对条件期望、条件方差、Poisson过程和Markov链的基本知识进行简单介绍之后,介绍了如何利用计算机产生随机数以及如何利用这些随机数产生任意分布的随机变量、随机过程等知识;介绍了一些分析统计数据的方法和技术...
  • 2.2.4 实际中如何用哈尔或沃尔什函数构建图像变换矩阵? 58 2.2.5 哈尔变换的基元图像看起来是什么样的? 61 2.2.6 可以定义元素仅为+1 或.1 的正交矩阵吗? 65 B2.5 对沃尔什函数的排列方式 65 2.2.7 哈达玛/...
  • 1.6 期望方差、协方差、相关系数 13 1.6.1 期望 13 1.6.2 方差 14 1.6.3 协方差 14 1.6.4 相关系数 15 第2章 机器学习基础 16 2.1 基本概念 16 2.1.1 大话机器学习本质 16 2.1.2 什么是神经网络 16 2.1.3 各种常见...
  • 统计学方法与数据分析(上下册)

    热门讨论 2013-12-29 11:32:47
    17.5计算期望均方的规则 17.6套抽样和裂区设计 17.7小结 补充练习 第十八章重复测量与交叉设计 18.1引言和案例 18.2有重复观测的单因子试验 18.3一个因子有重复观测的两因子试验 18.4交叉设计 18.5小结 ...
  • java 编程艺术

    2012-09-16 23:38:44
    9.4 计算达到某项期望年金所需的初始投资 278 9.5 根据给定投资计算年金的最大值 282 9.6 计算某项贷款的账户余额 286 9.7 创建用于金融类计算的Servlet 290 9.7.1 使用Tomcat 291 9.7.2 测试Servlet 291 ...
  • Java编程艺术 PDF

    2011-04-28 16:39:10
    9.4 计算达到某项期望年金所需的初始投资 278 9.5 根据给定投资计算年金的最大值 282 9.6 计算某项贷款的账户余额 286 9.7 创建用于金融类计算的Servlet 290 9.7.1 使用Tomcat 291 9.7.2 测试Servlet 291 9.7.3 把...
  • 1.0.7 如果一个传感器对应物理世界中的一个小片,如何能让多个传感器对应场景中的同一个小片?.................................................................2 1.0.8 什么是图像中一个像素位置亮度的物理含义...

空空如也

空空如也

1 2 3 4 5
收藏数 89
精华内容 35
关键字:

如何用期望计算方差