Марковские цепные стационарные распределения с scipy.sparse?

У меня цепь Маркова, данная как большая разреженная матрица A (Я построил матрицу в формате scipy.sparse.dok_matrix , но преобразование в другие или создание ее как csc_matrix прекрасное.)

Я хотел бы знать любое стационарное распределение p этой матрицы, которое является собственным вектором к собственному значению 1 . Все записи в этом собственном векторе должны быть положительными и содержать до 1, чтобы представить распределение вероятности.

Это означает, что я хочу любое решение для системы (AI) p = 0 , p.sum()=1 (где I=scipy.sparse.eye(*A.shape) – это матрица идентификаторов), но (AI) не будет иметь полный ранг, и даже вся система может быть недооценена. Кроме того, могут быть сгенерированы собственные векторы с отрицательными элементами, которые не могут быть нормализованы как допустимые распределения вероятностей. Предотвращение отрицательных записей в p было бы неплохо.

  • Использование scipy.sparse.linalg.eigen.eigs не является решением: оно не позволяет указывать аддитивное ограничение. (Нормализация не помогает, если в собственных векторах содержатся отрицательные записи.) Кроме того, он немного отличается от истинных результатов, а иногда имеет проблемы сходимости и ведет себя хуже, чем scipy.linalg.eig . (Кроме того, я использовал режим сдвигового инвертирования, который улучшает поиск типов собственных значений, которые я хочу, но не их качество. Если я его не использую, это еще больше перебор, потому что меня интересует только одно собственное значение. )

  • Преобразование в плотные матрицы и использование scipy.linalg.eig не является решением: помимо проблемы с отрицательной записью матрица слишком велика.

  • Использование scipy.sparse.spsolve не является очевидным решением: матрица либо не квадратная (при объединении аддитивного ограничения и условия собственного вектора), либо не полного ранга (при попытке указать их отдельно в некотором роде), иногда и ни один.

Есть ли хороший способ численно получить стационарное состояние цепи Маркова, заданное как разреженная матрица с использованием python? Если есть способ получить исчерпывающий список (и, возможно, также почти стационарные состояния), это оценено, но не обязательно.

  • Итерация через вектор scipy.sparse (или матрицу)
  • SciPy и scikit-learn - ValueError: несоответствие размеров
  • Многомерный нормальный CDF в Python с использованием scipy
  • Как извлечь точки из графика?
  • Более эффективно рисовать полигоны с помощью matplotlib
  • Раскрасить диаграмму Вороного
  • Scipy разреженные матрицы - назначение и использование различных реализаций
  • Проблема с типом данных с использованием scipy.spatial
  • 4 Solutions collect form web for “Марковские цепные стационарные распределения с scipy.sparse?”

    Несколько статей с резюме возможных подходов можно найти у Google ученого, вот что: http://www.ima.umn.edu/preprints/pp1992/932.pdf

    То, что сделано ниже, представляет собой комбинацию предложения @Helge Dietert выше о разделении на сильно связанные компоненты и подход № 4 в документе, приведенном выше.

     import numpy as np import time # NB. Scipy >= 0.14.0 probably required import scipy from scipy.sparse.linalg import gmres, spsolve from scipy.sparse import csgraph from scipy import sparse def markov_stationary_components(P, tol=1e-12): """ Split the chain first to connected components, and solve the stationary state for the smallest one """ n = P.shape[0] # 0. Drop zero edges P = P.tocsr() P.eliminate_zeros() # 1. Separate to connected components n_components, labels = csgraph.connected_components(P, directed=True, connection='strong') # The labels also contain decaying components that need to be skipped index_sets = [] for j in range(n_components): indices = np.flatnonzero(labels == j) other_indices = np.flatnonzero(labels != j) Px = P[indices,:][:,other_indices] if Px.max() == 0: index_sets.append(indices) n_components = len(index_sets) # 2. Pick the smallest one sizes = [indices.size for indices in index_sets] min_j = np.argmin(sizes) indices = index_sets[min_j] print("Solving for component {0}/{1} of size {2}".format(min_j, n_components, indices.size)) # 3. Solve stationary state for it p = np.zeros(n) if indices.size == 1: # Simple case p[indices] = 1 else: p[indices] = markov_stationary_one(P[indices,:][:,indices], tol=tol) return p def markov_stationary_one(P, tol=1e-12, direct=False): """ Solve stationary state of Markov chain by replacing the first equation by the normalization condition. """ if P.shape == (1, 1): return np.array([1.0]) n = P.shape[0] dP = P - sparse.eye(n) A = sparse.vstack([np.ones(n), dP.T[1:,:]]) rhs = np.zeros((n,)) rhs[0] = 1 if direct: # Requires that the solution is unique return spsolve(A, rhs) else: # GMRES does not care whether the solution is unique or not, it # will pick the first one it finds in the Krylov subspace p, info = gmres(A, rhs, tol=tol) if info != 0: raise RuntimeError("gmres didn't converge") return p def main(): # Random transition matrix (connected) n = 100000 np.random.seed(1234) P = sparse.rand(n, n, 1e-3) + sparse.eye(n) P = P + sparse.diags([1, 1], [-1, 1], shape=P.shape) # Disconnect several components P = P.tolil() P[:1000,1000:] = 0 P[1000:,:1000] = 0 P[10000:11000,:10000] = 0 P[10000:11000,11000:] = 0 P[:10000,10000:11000] = 0 P[11000:,10000:11000] = 0 # Normalize P = P.tocsr() P = P.multiply(sparse.csr_matrix(1/P.sum(1).A)) print("*** Case 1") doit(P) print("*** Case 2") P = sparse.csr_matrix(np.array([[1.0, 0.0, 0.0, 0.0], [0.5, 0.5, 0.0, 0.0], [0.0, 0.0, 0.5, 0.5], [0.0, 0.0, 0.5, 0.5]])) doit(P) def doit(P): assert isinstance(P, sparse.csr_matrix) assert np.isfinite(P.data).all() print("Construction finished!") def check_solution(method): print("\n\n-- {0}".format(method.__name__)) start = time.time() p = method(P) print("time: {0}".format(time.time() - start)) print("error: {0}".format(np.linalg.norm(PTdot(p) - p))) print("min(p)/max(p): {0}, {1}".format(p.min(), p.max())) print("sum(p): {0}".format(p.sum())) check_solution(markov_stationary_components) if __name__ == "__main__": main() 

    EDIT : обнаружена ошибка — csgraph.connected_components возвращает также чисто разлагающиеся компоненты, которые необходимо отфильтровать.

    Обращаясь к тому, что стационарные решения не уникальны и что решение не может быть неотрицательным.

    Это означает, что ваша цепь Маркова не является неприводимой и что вы можете разбить проблему на неприводимые цепи Маркова. Для этого вы хотите найти замкнутые коммуникационные классы вашей цепочки Маркова, которые по сути являются изучением связанных компонентов вашего переходного графа ( Wikipedia предлагает некоторые линейные алгоритмы для поиска сильно связанных компонентов). Более того, вы можете преодолеть все открытые коммуникационные классы, поскольку каждое стационарное состояние должно исчезнуть на них.

    Если у вас есть замкнутые коммуникационные классы C_1, …, C_n, ваша проблема, надеюсь, будет разбита на небольшие простые части: цепь Маркова на каждом замкнутом классе C_i теперь неприводима, так что ограниченная матрица перехода M_i имеет ровно один собственный вектор с собственным значением 0 и этот собственный вектор имеет только положительные компоненты (см. теорему Перрона-Фробениуса). Таким образом, мы имеем ровно одно стационарное состояние х_1.

    Стационарные состояния вашей полной цепи Маркова теперь представляют собой линейные комбинации x_i из ваших закрытых классов. На самом деле это все стационарные состояния.

    Чтобы найти стационарные состояния x_i, вы можете последовательно применять M_i, и итерация сходится к этому состоянию (это также сохранит вашу нормализацию). В общем, трудно сказать скорость конвергенции, но это дает вам простой способ повысить точность вашего решения и проверить решение.

    Это решение, возможно, underspecified Matrix уравнение, и поэтому может быть сделано с scipy.sparse.linalg.lsqr . Я не знаю, как убедиться, что все записи положительные, но, кроме того, это очень просто.

     import scipy.sparse.linalg states = A.shape[0] # I assume that the rows of A sum to 1. # Therefore, In order to use A as a left multiplication matrix, # the transposition is necessary. eigvalmat = (A - scipy.sparse.eye(states)).T probability_distribution_constraint = scipy.ones((1, states)) lhs = scipy.sparse.vstack( (eigvalmat, probability_distribution_constraint)) B = numpy.zeros(states+1) B[-1]=1 r = scipy.sparse.linalg.lsqr(lhs, B) # r also contains metadata about the approximation process p = r[0] 

    Используйте силовую итерацию (например): http://www.google.com/search?q=power%20iteration%20markov%20chain

    Или вы можете использовать режим переключения смены scipy.sparse.linalg.eig (который является ARPACK) для поиска собственных значений, близких к 1. «Указание» нормализации не требуется, так как вы можете просто нормализовать значения после этого.

    Python - лучший язык программирования в мире.