I would like to construct a square matrix which has the following form:
1.0, A, A^2, A^3, … A^n
0.0, 1.0, A, A^2, A^3, … A^n-1
0.0, 0.0, 1.0, A, A^2, … A^n-2
…
0.0, 0.0, 0.0, … 1.0
ie, the leading diagonal has ones, left/lower triangle is zeros, and right/upper triangle has power series, where ‘A’ is some random variable, eg pm.Normal()
the dimension of the matrix is determined by the input data.
You can use an elementwise exponentiation to do this. See if this code helps:
import pymc3 as pm
import theano.tensor as tt
import numpy as np
np.set_printoptions(precision=2)
# Dimension of square matrix
d = 5
C = np.zeros([d,d])
# Build a matrix in Numpy ahead of time
for i in range(d):
C[i, i:] = np.arange(1, (d-i)+1)
with pm.Model() as model:
A = pm.Normal('A')
matrix = pm.Deterministic('matrix',tt.triu(A ** pm.math.constant(C)))
trace = pm.sample()
trace['matrix'][0]