Python实现马尔可夫链入门指南

为什么要使用马尔可夫链？

马尔可夫链在数学中有着广泛的应用。它还广泛应用于经济学、博弈论、传播论、遗传学和金融学等领域。它经常出现在统计领域，尤其是贝叶斯统计和信息论领域。事实上，马尔可夫的连锁店为自动驾驶汽车控制系统、机场乘客排队、汇率等问题提供了解决方案。 PageRank，由第一个Google搜索引擎提出，是一种基于马尔可夫流程的算法。 Reddit 有一个名为模拟器 subreddit 的 subreddit。所有帖子和评论均由马尔可夫链自动创建。很不错！

马尔可夫链

马尔可夫链是一个具有马尔可夫属性的随机过程。随机过程或具有随机属性的数学对象是由一系列随机生成定义的。 马尔可夫 链有很多变体，因为它们具有离散状态空间（随机变量的一组可能值）或离散索引集（通常代表时间）。通常所说的“马尔可夫链”是指具有离散时间集合的过程，即单时间马尔可夫链（DTMC）。

Timeless Poutine Chain

特定时间马尔可夫链中包含的系统的所有阶段都处于一种状态，并且状态在阶段之间随机变化。这些步骤通常以时间来比较（但您也可以考虑物理距离或特定条件）。有限时间 Poutine 链是随机变量 X1, X2 的序列，概率的数学公式表示为：

Pr(Xn+1 = x | X1 = x1, X2 = x2, …, Xn = xn ) = Pr( Xn+1 = x | 概率+1只与之前的概率Xn有关。因此，需要知道之前的状态才能确定当前状态的概率分布，满足条件独立性（即：你需要知道当前状态才能确定下一个状态）。状态空间可以是任何东西：字母、数字、篮球得分。马尔可夫的状态链实现使用有限或不可数的状态空间，更容易进行统计分析。

模型

Poutine 链由概率自动化表示（相信我，它并不像听起来那么难！）。系统结构的变化称为过渡。每个状态发生变化的概率称为转移概率。概率自动化涉及将已知转换的概率转换为转换方程再转换为转换矩阵。

您还可以将 马尔可夫 链视为有向图，其中图的 n 条边由从时间 n 的状态移动到时间 n+1 的状态的概率来标记，Pr(Xn+1 = x | Xn = xn)。该公式可以理解为从已知状态Xn转变到状态Xn+1的概率。这个概念可以用每 n n+1 的 转移矩阵 来表示。状态空间中的每个状态首先作为转移矩阵的行出现，然后作为列出现。矩阵的每个元素代表从该行表示的状态转移到该列状态的概率。

如果马尔可夫链有N个状态，那么转移矩阵是N x N维，其中（I，J）表示从状态I到状态J的转移概率。另外，转移矩阵必须是概率矩阵，即每行各元素之和必须为1，为什么呢？因为每条线都代表了它自己的电位分布。

因此，模型的主要特征包括：状态空间、定义特定转移概率的转移矩阵以及初始分布提供的初始位置条件。

看起来很难？

让我们看一个简单的例子来帮助理解这些概念：

如果Cj很少刻薄，他会去跑步，吃很多冰淇淋（译者有福了：一开始的gooble应该是狼吞虎咽。），或者休息一下，休息一下调整一下。

根据过去的数据，如果他通过睡觉来调节心情，那么他第二天去跑步的可能性会增加60%，赖床的可能性会增加20%，吃大份的可能性会增加20%的食物。冰淇淋。

如果他去跑步休息，那么他第二天有60%的机会再跑步，有30%的机会他会吃冰淇淋，只有10%的机会他会睡觉。

最后，如果他在悲伤时沉迷于冰淇淋，那么只有10%的人第二天会继续吃冰淇淋，70%的人会跑步，20%的人会睡觉。

上面状态图中所示的马尔可夫链有3种可能的状态：睡眠、运行和冰冻。所以转移矩阵是一个3 x 3的矩阵。注意离开某个状态的箭头的值之和必须为1。这等于每行元素的和。状态矩阵的为​​1，表示概率分布。转移矩阵中每个元素的含义对于状态图中的每个状态都是相同的。

这个例子应该可以帮助您理解与马尔可夫链相关的一些不同概念。但这个理论如何应用于现实世界呢？

通过这个例子，你应该能够回答这类问题：“从睡眠状态开始，2天后Cj选择run（运行状态）的概率是多少？””

我们一起算算吧。要从睡眠状态转到运行状态，Cj 有以下选择：第一天继续睡眠，第二天运行（0.2 ⋅ 0.6）；第一天恢复跑步并在第二天恢复跑步 (0.6 ⋅ 0.6)； 1. 一天吃冰淇淋，第二天跑步（0.2 ⋅ 0.7）。计算出的概率为：((0.2 ⋅ 0.6) + (0.6 ⋅ 0.6) + (0.2 ⋅ 0.7)) = 0.62。所以，从睡眠状态开始，2天后Cj有62%的机会进入运行状态。

希望这个例子可以告诉你马尔可夫链网络可以解决哪些问题。

同时可以更好地理解马尔可夫链的几个重要属性：

• 互操作性：如果一个马尔可夫链可以从任何状态转移到任何状态，那么它就不能重复。换句话说，如果两个状态之间的步骤序列的概率为正，则它是不可逆的。
• 条件：如果马尔可夫链在大于1的整数倍后返回到相同状态，则马尔可夫链的状态是连续的。因此，从状态“i”开始，只有经过整数个周期“k”后才能返回到“i”，其中k是满足条件的所有数字中的最大值。如果 k = 1，则状态“i”不是周期性的，但如果 k > 1，则“i”是周期性的。
• 瞬态和递归性：如果从状态“i”开始并且可能无法返回到“i”，则“i”具有传递性。否则，它会重复（或永久）。如果一个状态可以在有限步内重复，则该状态是可重复的，否则是不可重复的。
• 遍历性：如果状态“i”满足非周期性和正递推，则它具有遍历性。如果不可约马尔可夫链中的每个状态都具有遍历性，那么这条马尔可夫链也具有遍历性。
• 避免状态：如果“i”无法移动到另一个状态，则“i”处于状态。因此，如果i≠j，pii = 1且pij = 0，则“i”处于被动状态。如果马尔可夫链中的所有状态都能达到拉动状态，则称为拉动状态的马尔可夫链。

提示：访问此网站以直观地了解马尔可夫的监管链。

使用Python实现马尔可夫链

让我们使用Python来实现上面的例子。当然，实际库实现的马尔可夫链的效率会高很多。这里有一些示例代码可以帮助您入门...

首先加载要使用的库。

import numpy as np
import random as rm
复制代码

然后定义状态及其可能性，即转移矩阵。请记住，由于存在三种状态，因此矩阵是 3 X 3 维。此外，还必须确定传输路径，这也可以在矩阵中显示。

# 状态空间
states = ["Sleep","Icecream","Run"]

# 可能的事件序列
transitionName = [["SS","SR","SI"],["RS","RR","RI"],["IS","IR","II"]]

# 概率矩阵（转移矩阵）
transitionMatrix = [[0.2,0.6,0.2],[0.1,0.6,0.3],[0.2,0.7,0.1]]
复制代码

别忘了，确保可能性的总和是 1。而且，写代码的时候多打印一些错误信息并没有什么问题！

if sum(transitionMatrix[0])+sum(transitionMatrix[1])+sum(transitionMatrix[1]) != 3:
    print("Somewhere, something went wrong. Transition matrix, perhaps?")
else: print("All is gonna be okay, you should move on!! ;)")
复制代码

All is gonna be okay, you should move on!! ;)
复制代码

现在让我们看看。我们将使用 numpy.random.choice 从可能的更改集中选择一个随机样本。代码中大部分参数的含义都可以通过参数名称找到，但参数p可能比较混乱。它是一个可选参数，可以传递给样本概率分布。这里引入了传递矩阵。

# 实现了可以预测状态的马尔可夫模型的函数。
def activity_forecast(days):
    # 选择初始状态
    activityToday = "Sleep"
    print("Start state: " + activityToday)
    # 应该记录选择的状态序列。这里现在只有初始状态。
    activityList = [activityToday]
    i = 0
    # 计算 activityList 的概率
    prob = 1
    while i != days:
        if activityToday == "Sleep":
            change = np.random.choice(transitionName[0],replace=True,p=transitionMatrix[0])
            if change == "SS":
                prob = prob * 0.2
                activityList.append("Sleep")
                pass
            elif change == "SR":
                prob = prob * 0.6
                activityToday = "Run"
                activityList.append("Run")
            else:
                prob = prob * 0.2
                activityToday = "Icecream"
                activityList.append("Icecream")
        elif activityToday == "Run":
            change = np.random.choice(transitionName[1],replace=True,p=transitionMatrix[1])
            if change == "RR":
                prob = prob * 0.5
                activityList.append("Run")
                pass
            elif change == "RS":
                prob = prob * 0.2
                activityToday = "Sleep"
                activityList.append("Sleep")
            else:
                prob = prob * 0.3
                activityToday = "Icecream"
                activityList.append("Icecream")
        elif activityToday == "Icecream":
            change = np.random.choice(transitionName[2],replace=True,p=transitionMatrix[2])
            if change == "II":
                prob = prob * 0.1
                activityList.append("Icecream")
                pass
            elif change == "IS":
                prob = prob * 0.2
                activityToday = "Sleep"
                activityList.append("Sleep")
            else:
                prob = prob * 0.7
                activityToday = "Run"
                activityList.append("Run")
        i += 1  
    print("Possible states: " + str(activityList))
    print("End state after "+ str(days) + " days: " + activityToday)
    print("Probability of the possible sequence of states: " + str(prob))

# 预测 2 天后的可能状态
activity_forecast(2)
复制代码

Start state: Sleep
Possible states: ['Sleep', 'Sleep', 'Run']
End state after 2 days: Run
Probability of the possible sequence of states: 0.12
复制代码

可以从睡眠状态及其潜力获得结果。进一步扩展此函数，您可以在睡眠状态下启动它并重复数百次以获得特定结束状态的概率。我们重写函数activity_forecast，添加一个循环...

def activity_forecast(days):
    # 选择初始状态
    activityToday = "Sleep"
    activityList = [activityToday]
    i = 0
    prob = 1
    while i != days:
        if activityToday == "Sleep":
            change = np.random.choice(transitionName[0],replace=True,p=transitionMatrix[0])
            if change == "SS":
                prob = prob * 0.2
                activityList.append("Sleep")
                pass
            elif change == "SR":
                prob = prob * 0.6
                activityToday = "Run"
                activityList.append("Run")
            else:
                prob = prob * 0.2
                activityToday = "Icecream"
                activityList.append("Icecream")
        elif activityToday == "Run":
            change = np.random.choice(transitionName[1],replace=True,p=transitionMatrix[1])
            if change == "RR":
                prob = prob * 0.5
                activityList.append("Run")
                pass
            elif change == "RS":
                prob = prob * 0.2
                activityToday = "Sleep"
                activityList.append("Sleep")
            else:
                prob = prob * 0.3
                activityToday = "Icecream"
                activityList.append("Icecream")
        elif activityToday == "Icecream":
            change = np.random.choice(transitionName[2],replace=True,p=transitionMatrix[2])
            if change == "II":
                prob = prob * 0.1
                activityList.append("Icecream")
                pass
            elif change == "IS":
                prob = prob * 0.2
                activityToday = "Sleep"
                activityList.append("Sleep")
            else:
                prob = prob * 0.7
                activityToday = "Run"
                activityList.append("Run")
        i += 1    
    return activityList

# 记录每次的 activityList
list_activity = []
count = 0

# `range` 从第一个参数开始数起，一直到第二个参数（不包含）
for iterations in range(1,10000):
        list_activity.append(activity_forecast(2))

# 查看记录到的所有 `activityList`    
#print(list_activity)

# 遍历列表，得到所有最终状态是跑步的 activityList
for smaller_list in list_activity:
    if(smaller_list[2] == "Run"):
        count += 1

# 计算从睡觉状态开始到跑步状态结束的概率
percentage = (count/10000) * 100
print("The probability of starting at state:'Sleep' and ending at state:'Run'= " + str(percentage) + "%")
复制代码

The probability of starting at state:'Sleep' and ending at state:'Run'= 62.419999999999995%
复制代码

那么问题来了，为什么计算结果会趋向于62%呢？

注意这是真正的“大数定律”在起作用。大数定律是概率论的一条定律，指出当试验次数较多时，相似事件的频率往往会发生变化。也就是说，随着测试次数的增加，实际率趋于理论或预测概率。

马尔可夫的想法

马尔可夫链教程就到此为止。本文介绍了马尔可夫及其资产的链条。 Simple 马尔可夫 Chains 是开始学习 Python 数据科学的安全方法。

作者：cdpath
来源：掘金