【发布时间】:2016-04-17 02:35:52
【问题描述】:
跟进this question 关于如何找到马尔可夫稳态,我现在遇到的问题是它在我的实验室计算机上完美运行,但在任何其他计算机上都无法运行。具体来说,它总是找到正确数量的近一特征值,因此哪些节点是吸引子节点,但它并不能始终找到所有这些节点,并且它们没有被正确分组。例如,使用下面的 64x64 转换矩阵,它不工作的计算机总是随机产生三个不同的错误吸引子集合之一。在下面较小的矩阵 M1 上,所有测试的计算机都得到相同的、正确的吸引子组和平稳分布的结果。
所有测试的机器都运行 Win7x64 和 WinPython-64bit-2.7.9.4。一台计算机总是做对,另外三台总是以同样的方式出错。根据我发现的几篇帖子like this 和this,这听起来可能是由于计算的浮点精度不同造成的。不幸的是,我不知道如何解决这个问题;我的意思是,我不知道如何更改从矩阵中提取左特征值的代码,以强制实现所有计算机都可以处理的特定精度(而且我认为它不必为此非常准确) .
这只是我目前对结果可能有何不同的最佳猜测。如果您对为什么会发生这种情况以及如何阻止这种情况发生有更好的了解,那就太好了。
如果有办法让 scipy 从运行到运行和计算机到计算机保持一致,我认为这不会取决于我的方法的细节,但因为它是被请求的,所以就在这里。两个矩阵都有 3 个吸引子。在 M1 中,第一个 [1,2] 是两个状态的轨道,另外两个 [7] 和 [8] 是平衡的。 M2 是一个64x64 transition matrix,在 [2] 和 [26] 处具有平衡,以及使用 [7,8] 的轨道。
但它有时会报告 [[26],[2],[26]] 有时会报告 [[2,7,8,26],[2],[26]] 而不是找到那组吸引子有时...每次运行都不会得到相同的答案,而且永远不会得到 [[2],[7,8],[26]] (以任何顺序)。
import numpy as np
import scipy.linalg
M1 = np.array([[0.2, 0.8, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
[1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
[0.6, 0.4, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
[0.0, 0.0, 0.2, 0.0, 0.1, 0.1, 0.3, 0.3],
[0.0, 0.0, 0.2, 0.2, 0.2, 0.2, 0.1, 0.1],
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.5],
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0],
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0]])
M2 = np.genfromtxt('transitionMatrix1.csv', delimiter=',')
# For easy switching
M = M2
# Confirm the matrix is a valid Markov transition matrix
#print np.sum(M,axis=1)
其余的代码与上一个问题中的代码相同,为方便起见,将其包含在此处。
#create a list of the left eigenvalues and a separate array of the left eigenvectors
theEigenvalues, leftEigenvectors = scipy.linalg.eig(M, right=False, left=True)
# for stationary distribution the eigenvalues and vectors are always real, and this speeds it up a bit
theEigenvalues = theEigenvalues.real
#print theEigenvalues
leftEigenvectors = leftEigenvectors.real
#print leftEigenvectors
# set how close to zero is acceptable as being zero...1e-15 was too low to find one of the actual eigenvalues
tolerance = 1e-10
# create a filter to collect the eigenvalues that are near enough to zero
mask = abs(theEigenvalues - 1) < tolerance
# apply that filter
theEigenvalues = theEigenvalues[mask]
# filter out the eigenvectors with non-zero eigenvalues
leftEigenvectors = leftEigenvectors[:, mask]
# convert all the tiny and negative values to zero to isolate the actual stationary distributions
leftEigenvectors[leftEigenvectors < tolerance] = 0
# normalize each distribution by the sum of the eigenvector columns
attractorDistributions = leftEigenvectors / leftEigenvectors.sum(axis=0, keepdims=True)
# this checks that the vectors are actually the left eigenvectors
attractorDistributions = np.dot(attractorDistributions.T, M).T
# convert the column vectors into row vectors (lists) for each attractor
attractorDistributions = attractorDistributions.T
print attractorDistributions
# a list of the states in any attractor with the stationary distribution within THAT attractor
#theSteadyStates = np.sum(attractorDistributions, axis=1)
#print theSteadyStates
【问题讨论】:
-
我给你留言了关于你在另一个问题上的回答,但是你有一个完整的复制案例可以添加吗?在最好的情况下,微妙的数值问题可能很难分析,但如果没有代码,就不可能提供任何有用的答案。
-
我支持@talonmies:请提供minimal reproducible example!
-
快速猜测:特征值和特征值没有排序返回。它们返回的顺序比值本身的变化要大得多。您是否有可能假设索引第一个/最后一个特征值/向量始终是最大或最小的?
-
当您说“特征向量错误”时,您是否检查过它们实际上不是特征向量? IE。
A x不等于lambda x在一定精度范围内吗? -
所有大写字母都没有帮助,你知道 :-)。我要问的是你是否检查过你没有得到退化特征值的特征向量的不同线性组合。
标签: python numpy scipy eigenvector eigenvalue