Generating random transition probability matrix python

Question

I was wondering whether there is an easy way to generate a square matrix filled with random numbers in python, given some conditions:

• the rows have to sum up to 1.
• the values on the diagonal should be significantly higher than the other values.

This is a stochastic matrix, and generating one is possible, however the tricky part is the condition about the diagonal. e.g. for a 4x4 matrix the output should look something like this:

[[0.90, 0.03, 0.03, 0.04],
[0.01, 0.98, 0.005, 0.005],
[0.04, 0.01, 0.92, 0.03],
[0.00, 0.02, 0.03, 0.95]]

Are there good ways to generate such matrix for a variable size?

Answer 1

Here is a quick-and-dirty solution

import random

k = 4

result = [[random.uniform(0, 0.1 / k) for i in range(k)] for j in range(k)]
for j, r in enumerate(result):
    r[j] += 1 - sum(r)

You may want to consider using a different random distribution and also have a look at numpy.

Answer 2

Here's a method with numpy.identity, starting with a k x k identity matrix, adding a drift term to it, and then normalizing.

import numpy as np

k = 4
result = np.identity(4)

# Add a random drift term.  We can guarantee that the diagonal terms
#     will be larger by specifying a `high` parameter that is < 1.
# How much larger depends on that term.  Here, it is 0.25.
result = result + np.random.uniform(low=0., high=.25, size=(k, k))

# Lastly, divide by row-wise sum to normalize to 1.
result = result / result.sum(axis=1, keepdims=1)

# Check
print(result)
print(result.sum(axis=1))

# [[ 0.80736896  0.00663004  0.06474194  0.12125906]
#  [ 0.03545472  0.79746194  0.10495657  0.06212678]
#  [ 0.08566011  0.02632533  0.79709851  0.09091605]
#  [ 0.07298408  0.05698381  0.1585878   0.71144431]]
#
# [ 1.  1.  1.  1.]

The above condensed to two lines:

result = np.identity(k) + np.random.uniform(low=0., high=.25, size=(k, k))
result /= result.sum(axis=1, keepdims=1)

Specifying a larger high parameter will get you a smaller ratio of "diagonals to the rest":

result = np.identity(k) + np.random.uniform(low=0., high=.60, size=(k, k))
result /= result.sum(axis=1, keepdims=1)
print(result.round(2))
# [[ 0.53  0.02  0.25  0.2 ]
#  [ 0.05  0.58  0.19  0.18]
#  [ 0.02  0.04  0.72  0.22]
#  [ 0.07  0.23  0.08  0.62]]