No way to extract the transition matrix from a Markov chain?

Issue #44 resolved
Paul Moore
created an issue

I'm experimenting with using Lea to analyze Markov chains. As an example, rules for rolling dice in certain games can be modelled as a Markov chain (a simplified example - roll 2 dice, rerolling doubles, is an absorbing chain where (n, m) is an absorbing state if n != m, and transitions to a new rolled pair if n == m).

One thing I'd like to do is extract the transition matrix from a `Chain` object. I can do this using the private `_next_state_lea_per_state` attribute, and deconstructing the `Alea` object for each starting state - but that's a little messy, and uses private information.

Is there a better way that I've missed? And if not, could one be added?

(Actually, Lea in general seems to make it hard to extract the underlying representation of an object - the best way I've found of getting the probability values for an `Alea` object is `list(zip(obj.support, obj.ps))` which is a little non-obvious (neither "support" nor "ps" are attribute names I would have immediately thought of).

1. repo owner

Hi Paul,

Happy to see your interest in Markov chain, a topic, which is for sure not much elaborated in Lea!

Let me first answer on the last point. Your need is addressed by methods described in the third part of the tutorial. You could for example replace your non-obvious expression by `obj.pmf_tuple` (well, you get a tuple, not a list!). I think that I should move this section in the first part of the wiki because it's seemingly a basic need of many users.

Now, about the main question, indeed, you didn’t miss anything. There is nothing out-of-the-box to get the transition matrix easily. Python is hopefully tolerant on accessing private attribute : as you have discovered already, you can get the transition NxN matrix as N tuple of N probabilities by

```tuple(a.ps for (_,a) in your_markov._next_state_lea_per_state)
```

and get the associated state names with

```tuple(s for (s,_) in your_markov._next_state_lea_per_state)
```

Definitely, I’m willing to add such access methods in Lea package. Now, is it exactly what you want? What method names would be meaningful? BTW, would it be sensible to rename `_next_state_lea_per_state` and promote it as public (and duly document it)?

Yes,… It’s now time to give your wish list. :)

2. reporter

I'm only really playing with Markov chains. There are some specific things I'm interested in, but they may not be particularly appropriate for Lea (and once I have the transition matrix, they aren't that hard to produce manually).

Specifically:

1. Is a Markov chain absorbing? (The ones I'm working with are, by construction, but checking that fact before I continue is not a bad idea).
2. Extract the number of transient states and the number of absorbing states.
3. Put the transition matrix into canonical form (trying to lay this out messes up Bitbucket's formatting, but basically 4 quadrants, Q in the upper left, R in the upper right, a zero matrix in the lower left and an identity matrix in the lower right - you probably already know this) and extract Q and R.
4. Calculate the fundamental matrix `(I-Q)**(-1)`

The matrix arithmetic really needs numpy (at least, it seems silly to avoid numpy) so would probably have to be optional in Lea anyway. So (4) probably isn't suitable for Lea.

Suggestions:

```# The name here is problematic as currently the states method does
# something different. But I can't think of a better name for a list of
# the states (as Python objects, not probability distributions)
all_states = chain.states()
# Boolean
chain.is_absorbing()
# A subset of all_states
absorbing = chain.absorbing_states()
transition_matrix = chain.transition_matrix()
# This is a bit long, not sure if it can be shortened?
q, r = chain.transition_matrix_canonical_decomposition()
```
3. repo owner

OK. That sounds interesting and general enough to be integrated in Lea. Let me work on this. Hopefully, I'll commit soon something that you could review.

4. repo owner

Document and change names of attributes/arguments of markov.Chain to be clearer (refs #44)

→ <<cset 4dd78701aecb>>

5. repo owner

Hi Paul,

Sorry for the long delay! I’ve eventually done something on the topic of absorbing MC (I know that you were able to manage, by hard ways, to get the information you need).

Everything is currently on the branch dev_upgrade_markov_chains of lea (to be merged on the main branch). There are 3 new methods on `markov.Chain` : `absorbing_mc_info`, `matrix` and `reachable_states`. The first is the main one, the last two are auxiliary but these could be helpful as such. These methods are documented in the code. Here is an example of use, for an absorbing MC:

```mc = markov.chain_from_matrix(
( 'A'  , 'B', 'C' ),
( 'A', ( 0.5 , 0.3 , 0.2 )),
( 'B', ( 0.0 , 1.0 , 0.0 )),
( 'C', ( 0.3 , 0.5 , 0.2 )))
print(mc.states)
# -> ('A', 'B', 'C')
print(mc.matrix())
# -> ((0.5, 0.3, 0.2), (0.0, 1.0, 0.0), (0.3, 0.5, 0.2))
print(mc.absorbing_mc_info())
# -> (True, ('A', 'C'), ('B',), ((0.5, 0.2), (0.3, 0.2)), ((0.3,), (0.5,)))
```

Notes:

1. Unlike your suggestions, I’ve preferred to return all data in one single method (`absorbing_mc_info`) because these are closely interdependent.
2. So far, numpy is not used in the methods; also, the fundamental matrix is not calculated. However, as I’ve done for other dependencies, I may envision an optional support to numpy, for example by adding an optional Boolean arg on `absorbing_mc_info` for asking to get Q/R matrices as numpy arrays and provide (I-Q)**(-1) at the end of the output tuple.
3. Of course, there remain work to do, at least: adding tests and documenting the new features on the wiki.

I would appreciate if you can give some feedback on these points before I go on.

6. repo owner

In markov.Chain, fix reachable_states method to remove duplicates and fix absorbing_mc_info to keep order of states (refs #44)

→ <<cset 4aa43be14e96>>

7. reporter

Hi, The above looks good to me. I haven't done much more since my original experiments - even though I had to play with internal attributes, that was enough to solve the problem I was working on at the time. But my solution would have been much cleaner with the new methods you're adding here. Optional numpy-based calculation of the fundamental matrix sounds like a reasonable way to go (although of course, it's possible to do that myself using the data you currently return in `absorbing_mc_info`).

8. repo owner

Thanks for your feedback, Paul. I continue to elaborate this way then.

9. repo owner

In markov.Chain, add as_array argument in matrix and absorbing_mc_info methods to return numpy arrays and calculate the fundamental matrix (refs #44)

→ <<cset 1b5fae309d7b>>

10. reporter

Looks awesome, thanks so much for doing this :-)