Defining a Custom PyMC Distribution with a Stochastic Matrix Variable

Question

I am trying to define the following distribution:

P(t) = exp(R*t)

Where R is 2x2 rates matrix that I want to solve for with the data (the sum of each of its rows must be 0). This is the log likelihood function, and some data:

import numpy as np
# The data is in the format: (row index, column index, t)
data = {(0,1,10), (1,0,20), (0,0,49), (1,1,12)}

def logp(values, ratesMatrix):
    r = []
    t = []
    for i in range(len(data)):
        r.append(ratesMatrix[values[i][0], values[i][1]])
        t.append(values[i][2])
    r = np.array(r, dtype=np.float64)
    t = np.array(t, dtype=np.float64)

    # log( prod( exp(r*t) ) )
    # == sum( log( exp(r*t) ) )
    # == sum( r*t )
    return np.sum(r*t)

How do I define the ratesMatrix variable as a stochastic 2x2 matrix where, each entry in the diagonal is between -1 and 0 (equally likely), and outside the diagonal between 0 and 1 (equally likely), such that the sum of each row is 0?

David David · Accepted Answer · 2018-08-19T22:20:11

I think this code works and is very simple:

import numpy as np

diagonal = np.random.randint(low=-1, high=1, size=2)
print(diagonal)

matrix = np.diag(diagonal)

if matrix[0, 0] == -1:
    matrix[0, 1] = 1
if matrix[1, 1] == -1:
    matrix[1, 0] = 1

print(matrix)

By the way, is that what you are looking for?

Defining a Custom PyMC Distribution with a Stochastic Matrix Variable

1 Answers