Generating random transition probability matrix python

3
votes

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?

2
pythonmatrix
How significantly higher? Do you have a probability distribution to quantify that? - Reti43
If all the values on the diagonal would be higher than 0.90, I would be happy! No probability distribution for quantifying yet. - DutchJ
Weeell, you could do random.uniform(0.9, 1.0) for the diagonal value and then generate a list of 3 other numbers and simply normalise them so that they sum up to 1 - diagonal_value. But I'm certain some statistician will point out how this method generates a biased sample of numbers, so you probably don't want to do it. By the way, if you're okay with using numpy, consider adding the tag. - Reti43
@DutchJ I can't imagine a situation where I would want random numbers but not even care what the distribution is. Can you give more context for the problem? - Max
@Reti43 the other way round may be ok. Generate all off-diagonal entries first using a distribution with an appropriately small mean and then fill in the diagonal - Paul Panzer

2 Answers

2
votes

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.

2
votes

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]]