Finde Markov Steady State mit linken Eigenwerten (mit numpy oder scipy)

Ich muss den stationären Zustand der Markov-Modelle mit den linken Eigenvektoren ihrer Übergangsmatrizen unter Verwendung eines Python-Codes finden.

Es ist bereits in dieser Frage festgestellt worden , dass scipy.linalg.eig nicht die tatsächlichen linken Eigenvektoren wie beschrieben liefert, sondern ein Fix wird dort demonstriert. Die offizielle Dokumentation ist meist nutzlos und unverständlich wie üblich.

Ein größeres Problem als das falsche Format ist, dass die erzeugten Eigenwerte nicht in einer bestimmten Reihenfolge (nicht sortiert und jeweils unterschiedlich) sind. Wenn du also die linken Eigenvektoren finden willst, die den 1 Eigenwerten entsprechen, musst du für sie jagen, und das stellt sein eigenes Problem dar (siehe unten). Die Mathematik ist klar, aber wie man Python bekommt, um dies zu berechnen und die richtigen Eigenvektoren zurückzugeben, ist nicht klar. Andere Antworten auf diese Frage, wie diese , scheinen nicht, die linken Eigenvektoren zu verwenden, also können die nicht korrekte Lösungen sein.

Diese Frage stellt eine Teillösung dar, aber sie berücksichtigt nicht die ungeordneten Eigenwerte größerer Übergangsmatrizen. Also, nur mit

leftEigenvector = scipy.linalg.eig(A,left=True,right=False)[1][:,0] leftEigenvector = leftEigenvector / sum(leftEigenvector) 

Ist nah, funktioniert aber nicht generell, weil der Eintrag in der Position [:,0] nicht der Eigenvektor für den korrekten Eigenwert sein kann (und in meinem Fall ist es normalerweise nicht).

Okay, aber die Ausgabe von scipy.linalg.eig(A,left=True,right=False) ist ein Array, in dem das [0] -Element eine Liste jedes Eigenwertes (nicht in beliebiger Reihenfolge) ist und in der Position gefolgt ist [1] durch ein Array von Eigenvektoren in einer entsprechenden Reihenfolge zu diesen Eigenwerten.

Ich kenne keinen guten Weg, diese ganze Sache durch die Eigenwerte zu sortieren oder zu suchen, um die richtigen Eigenvektoren herauszuziehen (alle Eigenvektoren mit Eigenwert 1, die durch die Summe der Vektoreinträge normalisiert sind). Mein Denken ist es, die Indizes der Eigenwerte zu erhalten Das gleich 1, und ziehe dann diese Spalten aus dem Array von Eigenvektoren. Meine Version davon ist langsam und umständlich. Zuerst habe ich eine Funktion (das funktioniert nicht ganz), um Positionen in einem letzten zu finden, der einem Wert entspricht:

 # Find the positions of the element a in theList def findPositions(theList, a): return [i for i, x in enumerate(theList) if x == a] 

Dann benutze ich es so, um die Eigenvektoren zu erhalten, die den Eigenwerten entsprechen = 1.

 M = transitionMatrix( G ) leftEigenvectors = scipy.linalg.eig(M,left=True,right=False) unitEigenvaluePositions = findPositions(leftEigenvectors[0], 1.000) steadyStateVectors = [] for i in unitEigenvaluePositions: thisEigenvector = leftEigenvectors[1][:,i] thisEigenvector / sum(thisEigenvector) steadyStateVectors.append(thisEigenvector) print steadyStateVectors 

Aber eigentlich funktioniert das nicht. Es gibt einen Eigenwert = 1.00000000e+00 +0.00000000e+00j , der nicht gefunden wird, obwohl zwei andere sind.

Meine Erwartung ist, dass ich nicht die erste Person bin, die Python benutzt, um stationäre Verteilungen von Markov-Modellen zu finden. Jemand, der mehr kompetent / erfahren ist, hat wahrscheinlich eine funktionierende allgemeine Lösung (ob mit numpy oder scipy oder nicht). Wenn man bedenkt, wie populäre Markov-Modelle ich erwartet habe, dort eine Bibliothek zu sein, nur für sie, um diese Aufgabe zu erfüllen, und vielleicht existiert es aber ich konnte keinen finden.

  • Wie erstelle ich eine Scipy Sparse Matrix aus einem Pandas Dataframe?
  • Nicht-NDFFrame-Objektfehler mit der Funktion pandas.SparseSeries.from_coo ()
  • Jakobische und hessische Eingaben in `scipy.optimize.minimize`
  • Informatikinformat in Python verarbeiten
  • Stop Scipy minimieren nach festgelegter Zeit
  • Speichern / laden scipy sparse csr_matrix im tragbaren Datenformat
  • Wie man elementweise auf einer Matrix des Typs scipy.sparse.csr_matrix operiert?
  • Scipy.optimize.fmin_cg: "'Gewünschter Fehler nicht unbedingt durch Präzisionsverlust erreicht'.
  • 2 Solutions collect form web for “Finde Markov Steady State mit linken Eigenwerten (mit numpy oder scipy)”

    Sie verknüpft Wie erfahre ich Eigenvektoren, die einem bestimmten Eigenwert einer Matrix entsprechen? Und sagte, es berechnet nicht den linken Eigenvektor, aber man kann das beheben, indem man mit der Transponierung arbeitet.

    Beispielsweise,

     In [901]: import numpy as np In [902]: import scipy.sparse.linalg as sla In [903]: M = np.array([[0.5, 0.25, 0.25, 0], [0, 0.1, 0.9, 0], [0.2, 0.7, 0, 0.1], [0.2, 0.3, 0, 0.5]]) In [904]: M Out[904]: array([[ 0.5 , 0.25, 0.25, 0. ], [ 0. , 0.1 , 0.9 , 0. ], [ 0.2 , 0.7 , 0. , 0.1 ], [ 0.2 , 0.3 , 0. , 0.5 ]]) In [905]: eval, evec = sla.eigs(MT, k=1, which='LM') In [906]: eval Out[906]: array([ 1.+0.j]) In [907]: evec Out[907]: array([[-0.32168797+0.j], [-0.65529032+0.j], [-0.67018328+0.j], [-0.13403666+0.j]]) In [908]: np.dot(evec.T, M).T Out[908]: array([[-0.32168797+0.j], [-0.65529032+0.j], [-0.67018328+0.j], [-0.13403666+0.j]]) 

    Um den Eigenvektor zu normalisieren (was du eigentlich real sein soll):

     In [913]: u = (evec/evec.sum()).real In [914]: u Out[914]: array([[ 0.18060201], [ 0.36789298], [ 0.37625418], [ 0.07525084]]) In [915]: np.dot(uT, M).T Out[915]: array([[ 0.18060201], [ 0.36789298], [ 0.37625418], [ 0.07525084]]) 

    Wenn Sie die Multiplizität des Eigenwertes 1 nicht scipy.linalg.eig , sehen Sie @ pv.'s Kommentar, der den Code mit scipy.linalg.eig . Hier ist ein Beispiel:

     In [984]: M Out[984]: array([[ 0.9 , 0.1 , 0. , 0. , 0. , 0. ], [ 0.3 , 0.7 , 0. , 0. , 0. , 0. ], [ 0. , 0. , 0.25, 0.75, 0. , 0. ], [ 0. , 0. , 0.5 , 0.5 , 0. , 0. ], [ 0. , 0. , 0. , 0. , 0. , 1. ], [ 0. , 0. , 0. , 0. , 1. , 0. ]]) In [985]: import scipy.linalg as la In [986]: evals, lvecs = la.eig(M, right=False, left=True) In [987]: tol = 1e-15 In [988]: mask = abs(evals - 1) < tol In [989]: evals = evals[mask] In [990]: evals Out[990]: array([ 1.+0.j, 1.+0.j, 1.+0.j]) In [991]: lvecs = lvecs[:, mask] In [992]: lvecs Out[992]: array([[ 0.9486833 , 0. , 0. ], [ 0.31622777, 0. , 0. ], [ 0. , -0.5547002 , 0. ], [ 0. , -0.83205029, 0. ], [ 0. , 0. , 0.70710678], [ 0. , 0. , 0.70710678]]) In [993]: u = lvecs/lvecs.sum(axis=0, keepdims=True) In [994]: u Out[994]: array([[ 0.75, -0. , 0. ], [ 0.25, -0. , 0. ], [ 0. , 0.4 , 0. ], [ 0. , 0.6 , 0. ], [ 0. , -0. , 0.5 ], [ 0. , -0. , 0.5 ]]) In [995]: np.dot(uT, M).T Out[995]: array([[ 0.75, 0. , 0. ], [ 0.25, 0. , 0. ], [ 0. , 0.4 , 0. ], [ 0. , 0.6 , 0. ], [ 0. , 0. , 0.5 ], [ 0. , 0. , 0.5 ]]) 

    Alles klar, ich musste bei der Umsetzung von Warrens Lösung einige Änderungen vornehmen. Es ist im Grunde das gleiche, so bekommt er den ganzen Kredit, aber die Realitäten der numerischen Approximationen mit numpy und scipy erforderten mehr Massage, die ich dachte, wäre hilfreich zu sehen, für andere versuchen, dies in der Zukunft zu tun. Ich habe auch die Variablen Namen zu super nervös gemacht.

    Bitte lassen Sie mich wissen, ob ich etwas falsches oder gibt es weitere empfohlene Verbesserungen (zB für Geschwindigkeit).

     # in this case my Markov model is a weighted directed graph, so convert that nx.graph (G) into it's transition matrix M = transitionMatrix( G ) #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 leftEigenvectors = leftEigenvectors.real # 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, but I guess it's not needed to usage #attractorDistributions = np.dot(attractorDistributions.T, M).T # convert the column vectors into row vectors (lists) for each attractor (the standard output for this kind of analysis) attractorDistributions = attractorDistributions.T # a list of the states in any attractor with the approximate stationary distribution within THAT attractor (eg for graph coloring) theSteadyStates = np.sum(attractorDistributions, axis=1) 

    Setzen Sie das alles in einem einfachen Copy-and-Paste-Format zusammen:

     M = transitionMatrix( G ) theEigenvalues, leftEigenvectors = scipy.linalg.eig(M, right=False, left=True) theEigenvalues = theEigenvalues.real leftEigenvectors = leftEigenvectors.real tolerance = 1e-10 mask = abs(theEigenvalues - 1) < tolerance theEigenvalues = theEigenvalues[mask] leftEigenvectors = leftEigenvectors[:, mask] leftEigenvectors[leftEigenvectors < tolerance] = 0 attractorDistributions = leftEigenvectors / leftEigenvectors.sum(axis=0, keepdims=True) attractorDistributions = attractorDistributions.T theSteadyStates = np.sum(attractorDistributions, axis=0) 

    Mit dieser Analyse auf einem generierten Markov-Modell produzierte ein Attraktor (von drei) mit einer stationären Verteilung von 0.19835218 und 0.80164782 im Vergleich zu den mathematisch genauen Werten von 0,2 und 0,8. Das ist also mehr als 0,1% Rabatt, eine Art großer Fehler für die Wissenschaft. Das ist kein echtes Problem, denn wenn Genauigkeit wichtig ist, dann, da die einzelnen Attraktoren identifiziert wurden, kann eine genauere Analyse des Verhaltens innerhalb jedes Attraktors unter Verwendung einer Matrix-Untermenge durchgeführt werden.

    Python ist die beste Programmiersprache der Welt.