【Python】【原创】马尔可夫链与动态行为和长期行为分析

白冥 · 发表于 2025-2-10 05:50:20

本帖最后由白冥于 2025-2-11 02:13 编辑

概述

      在生活中，一个离散系统可能有三种可能的结局，一种是永远在不停地变化，没有止息；一种是逐渐稳定，最终定格在遥远的明天；还有一种摇摆不定，振荡不止，似曾相识的昨日
不期而至。马尔可夫链便是这样的存在。
      本贴重点讨论马尔可夫链状态的可达性、类划分、首达概率、常返与瞬过、周期性等问题，并讨论平稳分布的概念。
      本帖代码经过重构，因此类结构和方法全部列出。代码和板书均由本人写就。

技术系列贴
      马尔可夫链（1）：马尔可夫链与马尔可夫性质

主要功能
      1)结构分解：使用互通性将转移图分解成多个互不干扰、彼此独立的环境进行局部性质的分析。
      2)动态行为分析：通过研究转移图的常返性和周期性间接测试系统的动态性质。
      3)长期稳态研究：通过局部规律模拟状态的长期分布。

类结构与方法
      Markov_chain 类
         ├ __init__ 方法
         ├ _build_P 方法
         ├ _get_classes 方法
         ├ is_irreducible 方法
         ├ _is_recurring_class 方法
         ├ get_recurring_classes 方法
         ├ _mul 方法
         ├ _get_period 方法
         ├ get_all_period_of_classes 方法
         └ get_stationary_distribution 方法

源代码

from math import gcd
from math import inf as INFTY
class Markov_chain:
def __init__(self, status:List[int], cpd:Dict[Tuple[int, int],float]):
self._status=status
self._P=self._build_P(status, cpd)
def _build_P(self, cpd):
n=len(status)
P=[[0]*n for _ in range(n)]
for e,p in cpd.items():
s_0,s_1=e
P[s_0][s_1]=p
return P
def _get_classes(self):
n = len(self._P)
reachable = [[self._P[i][j]>0 for j in range(n)] for i in range(n)]
for k in range(n):
for i in range(n):
for j in range(n):
reachable[i][j] = reachable[i][j] or (reachable[i][k] and reachable[k][j])
mutual = [[reachable[i][j] and reachable[j][i] for j in range(n)] for i in range(n)]
visited = [False]*n
classes = []
for i in range(n):
if visited[i]: continue
current_class=[i]
visited[i]=True
for j in range(n):
if mutual[i][j]:
visited[j]=True
current_class.append(j)
classes.append(current_class)
return classes
def is_irreducible(self):
return len(self._get_classes())==1
def _is_recurring_class(self, class_):
return not any(any(self._P[s_i][s_j]>0 for s_j in self._status if s_j not in class_)for s_i in class_)
def get_recurring_classes(self):
recurring_classes=[]
irreccurring_status=[]
for class_ in self._get_classes():
if self._is_recurring_class(class_):
recurring_classes.append(class_)
else:
for state in class_:
irreccurring_status.append(state)
return recurring_classes, irrecurring_status
def _mul(self, Q,P):
n=len(self._P)
return [[sum(Q[i][k]*P[k][j] for k in range(n))for j in range(n)] for i in range(n)]
def _get_period(self, class_):
s_0=class_[0]
n = len(self._P)
m = len(class_)
X=[]
P=[[p for p in r]for r in self._P]
Q=[[(1 if i=j else 0) for j in range(n)]for i in range(n)]
for k in range(1,m+1):
Q=self._mul(Q,P)
if Q[s_0][s_0]>0:
X.append(k)
return gcd(*X)
def get_all_period_of_classes(self):
recurring, irrecurring = self.get_recurring_classes()
period={}
for class_ in recurring:
period[tuple(class_)] = self._get_period(class_)
period[tuple(irrecurring)] = None
return period
def get_stationary_distribution(self,class_):
C=[[self._P[i][j] for j in class_] for i in class_]
n = len(C)
b = [0.0]*(n-1)+[1.0]
A=[]
for i in range(n-1):
A.append([-C[j][i] if i != j else 1 - C[j][i] for j in range(n)])
A.append([1.0]*n)
for i in range(n):
max_row = max(range(i, n), key=lambda r: abs(A[r][i]))
A[i], A[max_row] = A[max_row], A[i]
b[i], b[max_row] = b[max_row], b[i]
pivot = A[i][i]
for j in range(i, n):
A[i][j] /= pivot
b[i] /= pivot
for k in range(i+1, n):
factor = A[k][i]
for j in range(i, n):
A[k][j] -= factor * A[i][j]
b[k] -= factor * b[i]
pi = [0] * n
for i in range(n-1, -1, -1):
pi[i] = b[i]
for j in range(i+1, n):
pi[i] -= A[i][j] * pi[j]
return pi

复制代码

主要方法
1) _get_classes 方法
      功能：识别所有强连通类。
      步骤：这个方法的目标是将状态分成不同的强连通组件。首先，它初始化了一个reachable矩阵，记录每个状态i到j是否可达。然后通过三重循环进行传递闭包计算，这实际上是Floyd-Warshall算法的变种，用于检测可达性。之后，mutual矩阵用来记录i和j是否相互可达，从而确定强连通组件。然后通过遍历每个未被访问的节点，将同一类的节点收集起来。


2) _is_recurring_class方法
      功能：检查给定的类是否是常返的。
      步骤：这里判断的是类中的状态是否不会转移到类外的状态。具体来说，检查类中的任一状态s_i是否有任何转移到不在类中的s_j的概率大于0。如果有，则不是常返类。常返类意味着无法离开这个类。


3) _get_period 方法：
      功能：计算给定类的周期。
      步骤：这里选择类中的第一个状态s_0，然后计算从s_0出发返回s_0的所有可能步数的最大公约数。其中，只需要考虑最多m步（m是类的大小），因为超过m步的路径会有重复状态，从而可以分解成更小的环。


4) get_stationary_distribution 方法
      功能：用于计算平稳分布。
      步骤：首先将转移矩阵限制在类中的状态，然后构建增广矩阵A和b，用高斯消元法解平稳方程。


概念补充