在稀疏图中识别吸收 类
Identifying absorbing classes in a sparse graph
我的问题如下。我有几个吸收 类 的马尔可夫链,我想确定状态 space 中的哪些状态正在吸收,哪些是瞬态的。作为一个简单的例子,我们可以有以下稀疏转移矩阵,每行都有状态 0,1,2,3 的转移概率。
import numpy as np
from scipy.sparse import csr_matrix,csgraph
P = np.array([[0.6,0.4,0.,0.],[0.3,0.7,0.,0.],[0.2,0.4,0.2,0.2],[0.,0.,0.,1.]])
P = csr_matrix(P)
这里有两个吸收 类:状态 0-1 和状态 3。状态 2 是瞬态的。由于有 类 个吸收态,检查对角线上的概率 1 不足以识别吸收态 类。
我想知道如何以编程方式识别吸收 类。使用 csgraph
中的工具,我可以使用
找到强连通分量
components,labels = csgraph.connected_components(P,connection="strong")
print( components, labels )
3 [0 0 2 1]
所以,我有 3 个组件,每个组件都有自己的标签。吸收 类 肯定对应于这些组件之一,但也有瞬态组件。我如何找出哪个分量在吸收,哪个分量瞬态?
有人对此有好主意吗?
以下应该有效,注意一组吸收状态的总概率等于该组中的状态数:
components,labels = csgraph.connected_components(P, directed=True, connection='strong',return_labels=True)
absorbingStates = np.zeros(P.shape[0],dtype=bool)
for component in range(components):
indices = np.where(labels==component)[0]
n = len(indices)
if n==1:
probSum = P[indices,indices].sum() #np.ix_ gives an error if n==1
else:
probSum = P[np.ix_(indices,indices)].sum()
if np.isclose(probSum,n):
absorbingStates[indices] = True
我认为从 P
不断地 select 子矩阵不是很有效。欢迎更好的想法。
解决需要依赖浮点数学的马尔可夫链问题 (np.isclose
) 表明您的图形表示对于您尝试解决的问题不正确。
下面我介绍了一个使用 networkx 的替代解决方案(但我假设更改为 csgraph 并不难),它完全依赖于连接并计算出哪些边缘 link 强大的组件。
import networkx as nx
def absorbing_classes(g):
internal = set()
# Tarjan is approximately linear O(|V| + |E|)
cpts = list(nx.strongly_connected_component_subgraphs(g))
for sg in cpts:
for e in sg.edges():
internal.add(e)
# find all the edges that aren't part of the strongly connected components
# ~ O(E)
transient_edges = set(g.edges()) - internal
# find the start of the directed edge leading out from a component
# ~ O(E)
transient_srcs = set([ e[0] for e in transient_edges ])
# yield everything that don't have a vertex in transient_srcs
for sg in cpts:
if transient_srcs - set(sg.nodes()):
yield sg
例如,在您提供的图表上:
import numpy as np
P = np.array([[0.6,0.4,0.,0.],[0.3,0.7,0.,0.],[0.2,0.4,0.2,0.2],[0.,0.,0.,1.]])
g = nx.from_numpy_matrix(P, create_using=nx.DiGraph())
for subg in absorbing_classes(g):
print subg.nodes()
你得到两个吸收 类,其中一个节点 [0, 1]
,另一个节点 [3]
。
我的问题如下。我有几个吸收 类 的马尔可夫链,我想确定状态 space 中的哪些状态正在吸收,哪些是瞬态的。作为一个简单的例子,我们可以有以下稀疏转移矩阵,每行都有状态 0,1,2,3 的转移概率。
import numpy as np
from scipy.sparse import csr_matrix,csgraph
P = np.array([[0.6,0.4,0.,0.],[0.3,0.7,0.,0.],[0.2,0.4,0.2,0.2],[0.,0.,0.,1.]])
P = csr_matrix(P)
这里有两个吸收 类:状态 0-1 和状态 3。状态 2 是瞬态的。由于有 类 个吸收态,检查对角线上的概率 1 不足以识别吸收态 类。
我想知道如何以编程方式识别吸收 类。使用 csgraph
中的工具,我可以使用
components,labels = csgraph.connected_components(P,connection="strong")
print( components, labels )
3 [0 0 2 1]
所以,我有 3 个组件,每个组件都有自己的标签。吸收 类 肯定对应于这些组件之一,但也有瞬态组件。我如何找出哪个分量在吸收,哪个分量瞬态?
有人对此有好主意吗?
以下应该有效,注意一组吸收状态的总概率等于该组中的状态数:
components,labels = csgraph.connected_components(P, directed=True, connection='strong',return_labels=True)
absorbingStates = np.zeros(P.shape[0],dtype=bool)
for component in range(components):
indices = np.where(labels==component)[0]
n = len(indices)
if n==1:
probSum = P[indices,indices].sum() #np.ix_ gives an error if n==1
else:
probSum = P[np.ix_(indices,indices)].sum()
if np.isclose(probSum,n):
absorbingStates[indices] = True
我认为从 P
不断地 select 子矩阵不是很有效。欢迎更好的想法。
解决需要依赖浮点数学的马尔可夫链问题 (np.isclose
) 表明您的图形表示对于您尝试解决的问题不正确。
下面我介绍了一个使用 networkx 的替代解决方案(但我假设更改为 csgraph 并不难),它完全依赖于连接并计算出哪些边缘 link 强大的组件。
import networkx as nx
def absorbing_classes(g):
internal = set()
# Tarjan is approximately linear O(|V| + |E|)
cpts = list(nx.strongly_connected_component_subgraphs(g))
for sg in cpts:
for e in sg.edges():
internal.add(e)
# find all the edges that aren't part of the strongly connected components
# ~ O(E)
transient_edges = set(g.edges()) - internal
# find the start of the directed edge leading out from a component
# ~ O(E)
transient_srcs = set([ e[0] for e in transient_edges ])
# yield everything that don't have a vertex in transient_srcs
for sg in cpts:
if transient_srcs - set(sg.nodes()):
yield sg
例如,在您提供的图表上:
import numpy as np
P = np.array([[0.6,0.4,0.,0.],[0.3,0.7,0.,0.],[0.2,0.4,0.2,0.2],[0.,0.,0.,1.]])
g = nx.from_numpy_matrix(P, create_using=nx.DiGraph())
for subg in absorbing_classes(g):
print subg.nodes()
你得到两个吸收 类,其中一个节点 [0, 1]
,另一个节点 [3]
。