5. Dynamic Bayesian Network, Markov Chain

Let’s see how we can represent a Markov Chain (MC) as a Dynamic Bayesian Network (DBN). We will verify our results with the stationary distribution or steady state.

5.1. Model

The model, $$P$$, has two states, sunny and rainy. Our initial state, $$s$$, will be sunny.

:
import numpy as np

P = np.array([[0.9, 0.1], [0.5, 0.5]])
s = np.array([1.0, 0.0])

5.2. Power Method

We can use the Power Method to find the steady state.

$$s^n = s^0P^n$$

We will continue to compute $$s^n$$ until $$||P^n - P^{n-1}|| \lt \delta$$.

:
s.dot(P)
:
array([0.9, 0.1])
:
s.dot(P).dot(P)
:
array([0.86, 0.14])
:
import pandas as pd

P_curr = s.dot(P)
trace_df = [list(P_curr) + [np.nan]]

for it in range(100):
P_next = P_curr.dot(P)
d = np.linalg.norm(P_next - P_curr, 1)

if d < 0.0001:
break

trace_df.append(list(P_next) + [d])
P_curr = P_next

trace_df = pd.DataFrame(trace_df, columns=['sunny', 'rainy', 'delta'])
P_next
:
array([0.83335081, 0.16664919])

Here is the trace of $$s^n$$ until convergence.

:
trace_df
:
sunny rainy delta
0 0.900000 0.100000 NaN
1 0.860000 0.140000 0.080000
2 0.844000 0.156000 0.032000
3 0.837600 0.162400 0.012800
4 0.835040 0.164960 0.005120
5 0.834016 0.165984 0.002048
6 0.833606 0.166394 0.000819
7 0.833443 0.166557 0.000328
8 0.833377 0.166623 0.000131

5.3. Eigen decomposition

We can also find the steady state directly through eigen decomposition. The steady state is the eigenvector associated with the eigenvalue closest to 1.

:
def get_stationary_p(X):
S, U = np.linalg.eig(X.T)
r = (U[:,np.isclose(S, 1)][:,0] / U[:,np.isclose(S, 1)][:,0].sum()).real
return r

get_stationary_p(P)
:
array([0.83333333, 0.16666667])

5.4. Dynamic Bayesian Network

Another way to represent a MC is with a DBN. In general, the DBN structure consists of an endless number of nodes corresponding to time slices. However, in a DBN, we can unroll the DBN structure to just 2 consecutive time slices. It is assumed that the intra- (within a time slice) and inter-relationships (between two consecutive time slices) do not change over time.

:
from pybbn.graph.dag import Bbn
from pybbn.graph.edge import Edge, EdgeType
from pybbn.graph.jointree import EvidenceBuilder
from pybbn.graph.node import BbnNode
from pybbn.graph.variable import Variable
from pybbn.pptc.inferencecontroller import InferenceController

# create the nodes
w_curr = BbnNode(Variable(0, 'weather', ['rainy', 'sunny']), [1.0, 0.0])
w_next = BbnNode(Variable(2, 'weather_next', ['rainy', 'sunny']), [0.9, 0.1, 0.5, 0.5])

# create the network structure
bbn = Bbn() \

# convert the BBN to a join tree
join_tree = InferenceController.apply(bbn)

for node in [w_curr, w_next]:
potential = join_tree.get_bbn_potential(node)
print(node)
print(potential)
print('-' * 15)
0|weather|rainy,sunny
0=rainy|1.00000
0=sunny|0.00000
---------------
2|weather_next|rainy,sunny
2=rainy|0.90000
2=sunny|0.10000
---------------

We can perform inference on the DBN and use the probability vector of the next state as input to the following time slice.

:
import pandas as pd

def get_probs(jt, node):
entries = jt.get_bbn_potential(node).entries
return np.array([e.value for e in entries])

P_curr = get_probs(join_tree, w_curr)
jt = InferenceController.reapply(join_tree, {w_curr.id: list(P_curr)})

trace_df = [list(P_curr) + [np.nan]]

for _ in range(20):
P_next = get_probs(jt, w_next)
d = np.linalg.norm(P_next - P_curr, 1)

if d < 0.0001:
break

trace_df.append(list(P_next) + [d])
P_curr = P_next

jt = InferenceController.reapply(jt, {w_curr.id: list(P_curr)})

trace_df = pd.DataFrame(trace_df, columns=['sunny', 'rainy', 'delta'])
P_next
:
array([0.83335081, 0.16664919])
:
trace_df
:
sunny rainy delta
0 1.000000 0.000000 NaN
1 0.900000 0.100000 0.200000
2 0.860000 0.140000 0.080000
3 0.844000 0.156000 0.032000
4 0.837600 0.162400 0.012800
5 0.835040 0.164960 0.005120
6 0.834016 0.165984 0.002048
7 0.833606 0.166394 0.000819
8 0.833443 0.166557 0.000328
9 0.833377 0.166623 0.000131