Wie kann ich ein diskretes Markov-Modell mit Pymc erstellen?

Question

Ich versuche herauszufinden, wie man ein Markov-Kettenmodell mit diskreten Zuständen mit pymc richtig herstellt.

Als Beispiel (sehen Sie in nbviewer), machen wir eine Kette der Länge T = 10, wo der Markov-Zustand binär ist, die anfängliche Zustandsverteilung ist [0,2, 0,8] und dass die Wahrscheinlichkeit von Schaltzuständen in Zustand 1 ist 0,01, während im Zustand 2 ist 0,5

import numpy as np 
import pymc as pm 
T = 10 
prior0 = [0.2, 0.8] 
transMat = [[0.99, 0.01], [0.5, 0.5]] 

Um das Modell zu machen, mache ich eine Reihe von Zustandsvariablen, und eine Reihe von Übergangswahrscheinlichkeiten, die auf den Zustandsvariablen abhängen, (unter Verwendung der Funktion pymc.Index)

states = np.empty(T, dtype=object) 
states[0] = pm.Categorical('state_0', prior0) 
transPs = np.empty(T, dtype=object) 
transPs[0] = pm.Index('trans_0', transMat, states[0]) 

for i in range(1, T): 
    states[i] = pm.Categorical('state_%i' % i, transPs[i-1]) 
    transPs[i] = pm.Index('trans_%i' %i, transMat, states[i]) 

Probenahme das Modell zeigt, dass th e Staaten Rn sind, was sollten sie

model = pm.MCMC([states, transPs]) 
model.sample(10000, 5000) 
[np.mean(model.trace('state_%i' %i)[:]) for i in range(T)]  

druckt (mit Modell mit Kevin Murphy BNT Paket in Matlab gebaut verglichen) sein:

[-----------------100%-----------------] 10000 of 10000 complete in 7.5 sec 

[0.80020000000000002, 
0.39839999999999998, 
0.20319999999999999, 
0.1118, 
0.064199999999999993, 
0.044600000000000001, 
0.033000000000000002, 
0.026200000000000001, 
0.024199999999999999, 
0.023800000000000002] 

Meine Frage ist - das scheint nicht, wie die eleganteste Weg, eine Markov-Kette mit Pymc zu bauen. Gibt es einen saubereren Weg, der das Array von deterministischen Funktionen nicht benötigt?

Mein Ziel ist es, ein Pymc-Paket für allgemeinere dynamische Bayes-Netzwerke zu schreiben.

Answer 1

Soweit ich weiß, müssen Sie die Verteilung jedes Zeitschritts als deterministische Funktion des vorherigen Zeitschritts kodieren, denn das ist es - es gibt keine Zufälligkeit in den Parametern, weil Sie sie im Problem definiert haben Konfiguration. Ich denke jedoch, dass Ihre Frage vielleicht mehr darauf gerichtet war, eine intuitivere Darstellung des Modells zu finden. Ein alternativer Weg wäre, die Zeitschrittübergänge als eine Funktion des vorherigen Zeitschritts direkt zu codieren.

from pymc import Bernoulli, MCMC 

def generate_timesteps(N,p_init,p_trans): 
    timesteps=np.empty(N,dtype=object) 
    # A success denotes being in state 2, a failure being in state 1 
    timesteps[0]=Bernoulli('T0',p_init) 
    for i in xrange(1,N): 
     # probability of being in state 1 at time step `i` given time step `i-1` 
     p_i = p_trans[1]*timesteps[i-1]+p_trans[0]*(1-timesteps[i-1]) 
     timesteps[i] = Bernoulli('T%d'%i,p_i) 
    return timesteps 

timesteps = generate_timesteps(10,0.8,[0.001,0.5]) 
model = MCMC(timesteps) 
model.sample(10000) # no burn in necessary since we're sampling directly from the distribution 
[np.mean(model.trace(t).gettrace()) for t in timesteps] 

Answer 2

Falls Sie auf lange Sicht Verhalten Ihrer Markov-Kette aussehen wollen, kann das discreteMarkovChain Paket nützlich sein. Die Beispiele zeigen einige Ideen zum Aufbau einer diskreten Zustands-Markov-Kette durch Definieren einer Übergangsfunktion, die Ihnen für jeden Zustand die erreichbaren Zustände im nächsten Schritt und ihre Wahrscheinlichkeiten angibt. Sie könnten dieselbe Funktion verwenden, um den Prozess zu simulieren.