马尔可夫链在蒙特卡洛方法中的直观使用

容易吗? 我试过了


Netmology 数据科学讲师DomKlik的数据开发和工作主管Alexey Kuzmin翻译 Rahul Agarwal 的文章,内容涉及蒙特卡罗方法如何与马尔可夫链一起解决大型状态空间的问题。

与数据科学相关的每个人都曾经听说过使用马尔可夫链(MCMC)的蒙特卡洛方法。 在研究贝叶斯统计量时,有时会涉及到该主题,有时在使用诸如Prophet之类的工具时会涉及到该话题。

但是MCMC很难理解。 每次阅读这些方法时,我都注意到MCMC的本质隐藏在数学噪声的深层中,并且很难弄清这种噪声的背后。 我不得不花很多时间来理解这个概念。

在本文中,可以尝试用马尔可夫链来解释蒙特卡洛方法,以便清楚地了解它们的用途。 在下一篇文章中,我将重点介绍使用这些方法的更多方法。

因此,让我们开始吧。 MCMC由两个术语组成:蒙特卡洛链和马尔可夫链。 让我们谈谈每个。

蒙特卡洛




用最简单的术语来说,蒙特卡洛方法可以定义为简单的模拟。

蒙特卡洛方法(Monte Carlo Methods)的名字来自摩纳哥的蒙特卡洛赌场(Monte Carlo Casino)。 在许多纸牌游戏中,您需要知道赢得庄家的可能性。 有时,此概率的计算在数学上可能是复杂的或棘手的。 但是我们总是可以运行一次计算机模拟来多次玩整个游戏,并将概率视为胜利数除以游戏数。
这就是您需要了解的有关蒙特卡洛方法的全部信息。 是的,这只是一个简单的建模技术,名字很漂亮。

马尔可夫链




由于术语MCMC由两部分组成,因此您仍然需要了解什么是马尔可夫链 。 但是在进入马尔可夫链之前,让我们先谈谈马尔可夫特性。

假设有一个由M个可能的状态组成的系统,并且您从一种状态转移到另一种状态。 现在还不要让您感到困惑。 这种系统的一个特定示例是天气,天气从炎热到寒冷再到中度。 另一个例子是股票市场,它从看跌状态过渡到看涨停滞状态。

马尔可夫性质表明,对于在特定时间处于状态X n的给定过程,概率X n +1 = k(其中k是该过程可以进入的M个状态中的任何一个)仅取决于目前这种情况是什么。 而不是关于它如何达到当前状态。
用数学术语,我们可以用以下公式的形式表示:

为了清楚起见,您不必担心市场变得看涨的条件顺序。 下一个状态“看涨”的可能性仅取决于市场当前处于“看涨”状态的事实。 在实践中也很有意义。

具有马尔可夫特性的过程称为马尔可夫过程。 马尔可夫链为何重要? 由于其固定分布。

什么是固定分布?


在下面的示例中,我将尝试通过计算平稳分布来进行解释。 假设您有一个股市的马尔可夫过程,如下所示。

您具有一个转换概率矩阵,该矩阵确定从状态X i到X j转换的概率。

转移概率矩阵,Q

在瞬态概率矩阵Q中,给定当前状态“牛” = 0.9,下一个状态为“牛”的概率; 如果当前状态为“牛”,则下一个状态为“看涨”的概率为0.075。 依此类推。

好吧,让我们从某种特定状态开始。 我们的状态将由矢量[牛,熊,停滞]设置。 如果我们从“看跌”状态开始,向量将是这样的:[0,1,0]。 我们可以通过将当前状态向量乘以转移概率矩阵来计算下一个状态的概率分布。

请注意,概率加起来为1。

可以通过以下公式找到状态的以下分布:


依此类推。 最后,您将达到一个稳定的状态,其中该状态稳定下来:


对于上述的转移概率矩阵Q,平稳分布s为

您可以使用以下代码获得固定分布:

Q = np.matrix([[0.9,0.075,0.025],[0.15,0.8,0.05],[0.25,0.25,0.5]]) init_s = np.matrix([[0, 1 , 0]]) epsilon =1 while epsilon>10e-9:    next_s = np.dot(init_s,Q)    epsilon = np.sqrt(np.sum(np.square(next_s - init_s)))    init_s = next_s print(init_s) ------------------------------------------------------------------ matrix([[0.62499998, 0.31250002, 0.0625  ]]) 

您还可以从任何其他状态开始-实现相同的固定分布。 如果要确保这一点,请更改代码中的初始状态。

现在我们可以回答为什么固定分布如此重要的问题。

平稳分布很重要,因为它可以用来确定系统随机处于某个状态的概率。

对于我们的示例,可以说在62.5%的情况下,市场将处于“看涨”状态,31.25%处于“看跌”状态,而停滞则为6.25%。

直观地,您可以将其视为在链中随机游荡。


随机漫步

您在某个点上并选择下一个状态,同时考虑当前状态,观察下一个状态的概率分布。 根据这些节点的概率,我们可以比其他节点访问更多的节点。

这就是Google在互联网诞生之初解决搜索问题的方式。 问题是根据页面的重要性对页面进行排序。 Google使用Pagerank算法解决了该问题。 Google Pagerank算法应将状态视为页面,并将页面在固定分布中的概率视为其相对重要性。

现在我们直接考虑MCMC方法。

什么是带马尔可夫链的蒙特卡洛方法(MCMC)


在回答什么是MCMC之前,我先问一个问题。 我们了解Beta分布。 我们知道它的概率密度函数。 但是我们可以从这个分布中抽取样本吗? 您能提出一种方法吗?


想...

MCMC允许您从任何概率分布中进行选择。 当您需要从后验分布中进行选择时,这一点尤其重要。

该图显示了贝叶斯定理。

例如,您需要根据后验分布制作样本。 但是,计算后验分量和归一化常数(证据)是否容易? 在大多数情况下,您可以以可能性和先验概率乘积的形式找到它们。 但是计算归一化常数(p(D))不起作用。 怎么了 让我们仔细看看。

假设H仅采用3个值:

p(D)= p(H = H1).p(D | H = H1)+ p(H = H2).p(D | H = H2)+ p(H = H3).p(D | H = H3)

在这种情况下,p(D)易于计算。 但是,如果H的值是连续的,该怎么办? 可以很容易地计算出这个值,尤其是当H取无穷大时吗? 为此,必须解决一个复杂的积分。

我们要从后验分布中进行随机选择,但也要考虑p(D)为常数。

维基百科写道

带有马尔可夫链的蒙特卡洛方法是一类算法,用于根据马尔可夫链的构造从概率分布中进行采样,该马尔可夫链作为平稳分布具有所需的形式。 然后将一系列步骤后的链状状态用作所需分布的选择。 采样质量随步数的增加而提高。

让我们来看一个例子。 假设您需要Beta发行版中的样本。 其密度:


其中C是归一化常数。 实际上,这是α和β的某个函数,但是我想证明对于beta分布中的样本并不需要它,因此我们将其视为一个常数。

Beta分布问题确实很难解决,即使实际上不能解决。 实际上,您可能必须使用更复杂的分布函数,并且有时您将不知道标准化常数。

鉴于我们可以从均匀分布中进行选择(相对简单),MCMC方法通过提供可以创建具有β分布作为平稳分布的马尔可夫链的算法,使生活更轻松。

如果我们从一个随机状态开始,并根据某种算法几次进入下一个状态,则最终将创建一个具有β分布的Markov链作为平稳分布。 而且,我们长期以来发现的状态可以用作beta分布的样本。

这些MCMC算法之一是Metropolis-Hastings算法。

Metropolis-Hastings算法




直觉:


那么目的是什么?

直观地讲,我们希望沿着某个表面(我们的马尔可夫链)行走,这样我们在每个位置花费的时间与该位置的表面高度(要从中进行选择的期望概率密度)成正比。

因此,例如,我们希望在100米高的山顶上花费比在相邻50米山上花费更多的时间。 即使我们不知道曲面上点的绝对高度,也可以这样做:您只需要知道相对高度即可。 例如,如果山顶A的高度是山顶B的两倍,那么我们希望在A上花费的时间是在B上花费的两倍。

对于采用新的场所和规则提出了更为复杂的方案,但是主要思想如下:

  1. 选择一个新的“建议”位置。
  2. 找出此位置与当前位置相比是高还是低。
  3. 留在原地或搬到新的地方,其概率与这些地方的高度成比例。

MCMC的目的是从某种概率分布中进行选择,而不必知道其在任何点的确切高度(无需知道C)。
如果“漫游”过程设置正确,则可以确保达到此比例(花费的时间和分布高度之间)

算法:


现在,让我们以更正式的术语来定义和描述任务。 令s =(s1,s2,...,sM)为所需的平稳分布。 我们想创建一个具有这种平稳分布的马尔可夫链。 我们从带有过渡矩阵P的M状态的任意马尔可夫链开始,所以pij表示从状态i过渡到j的概率。

直观上,我们知道如何漫游马尔可夫链,但是马尔可夫链没有所需的平稳分布。 该链具有一些固定分布(我们不需要)。 我们的目标是改变我们在马尔可夫链上徘徊的方式,以使链具有所需的平稳分布。

为此:

  1. 从一个随机的初始状态开始。
  2. 通过查看转移矩阵P第i行中的转移概率,随机选择一个新的假定状态。
  3. 计算一种称为决策概率的度量,其定义为:aij = min(sj.pji / si.pij,1)。
  4. 现在扔一个硬币,它以aij的概率落在鹰的表面。 如果老鹰掉落,则接受要约,即进入下一个状态,否则,拒绝要约,即保持当前状态。
  5. 重复很多次。

经过大量测试,该链将收敛并具有固定分布s。 然后,我们可以将链状状态用作任何分布的样本。

通过执行此操作来采样beta分布,您唯一必须使用概率密度的方法是搜索决策的概率。 为此,将sj除以si(即归一化常数C被取消)。

Beta选择




现在我们来谈谈从beta分布采样的问题。

Beta分布是[0,1]上的连续分布,并且可以在[0,1]上具有无限值。 假设在[0,1]上具有无限状态的任意马尔可夫链P具有一个过渡矩阵P,使得pij = pji =矩阵中的所有元素。

我们不需要矩阵P,正如我们稍后将要看到的那样,但是我希望问题的描述尽可能与我们提出的算法相似。

  • 从从(0,1)上的均匀分布获得的随机初始状态i开始。
  • 通过查看转移矩阵P的第i行中的转移概率,随机选择一个新的假定状态。假定我们选择另一个状态Unif(0,1)作为假定状态j。
  • 计算测度,称为决策概率:


简化为:

由于pji = pij,所以

  • 现在扔硬币。 很有可能aij鹰会掉落。 如果老鹰掉落,那么您应该接受要约,即转移到下一个状态。 否则,值得拒绝报价,即保持相同状态。
  • 重复测试多次。

代码:


现在是时候从理论转向实践了。 我们将使用Python编写Beta版示例。

 impo rt rand om # Lets define our Beta Function to generate s for any particular state. We don't care for the normalizing constant here. def beta_s(w,a,b): return w**(a-1)*(1-w)**(b-1) # This Function returns True if the coin with probability P of heads comes heads when flipped. def random_coin(p): unif = random.uniform(0,1) if unif>=p: return False else: return True # This Function runs the MCMC chain for Beta Distribution. def beta_mcmc(N_hops,a,b): states = [] cur = random.uniform(0,1) for i in range(0,N_hops): states.append(cur) next = random.uniform(0,1) ap = min(beta_s(next,a,b)/beta_s(cur,a,b),1) # Calculate the acceptance probability if random_coin(ap): cur = next return states[-1000:] # Returns the last 100 states of the chain 

将结果与实际的beta分布进行比较。

 impo rt num py as np import pylab as pl import scipy.special as ss %matplotlib inline pl.rcParams['figure.figsize'] = (17.0, 4.0) # Actual Beta PDF. def beta(a, b, i): e1 = ss.gamma(a + b) e2 = ss.gamma(a) e3 = ss.gamma(b) e4 = i ** (a - 1) e5 = (1 - i) ** (b - 1) return (e1/(e2*e3)) * e4 * e5 # Create a function to plot Actual Beta PDF with the Beta Sampled from MCMC Chain. def plot_beta(a, b): Ly = [] Lx = [] i_list = np.mgrid[0:1:100j] for i in i_list: Lx.append(i) Ly.append(beta(a, b, i)) pl.plot(Lx, Ly, label="Real Distribution: a="+str(a)+", b="+str(b)) pl.hist(beta_mcmc(100000,a,b),normed=True,bins =25, histtype='step',label="Simulated_MCMC: a="+str(a)+", b="+str(b)) pl.legend() pl.show() plot_beta(0.1, 0.1) plot_beta(1, 1) plot_beta(2, 3) 




如您所见,这些值与beta分布非常相似。 因此,MCMC网络已达到稳定状态

在上面的代码中,我们创建了一个beta采样器,但是相同的概念适用于我们要从中进行选择的任何其他发行版。

结论




这是一个很棒的帖子。 祝贺您阅读到最后。

本质上,MCMC方法可能很复杂,但是它们为我们提供了极大的灵活性。 您可以通过MCMC使用选择从任何分配功能中选择。 通常,这些方法用于从后验分布中采样。

您还可以使用MCMC解决状态空间较大的问题。 例如,在背包问题或用于解密中。 在下一篇文章中,我将尝试为您提供更多有趣的示例。 请继续关注。

来自编辑


Source: https://habr.com/ru/post/zh-CN460497/


All Articles