# Learn linear transformation with pymc

Hello!
I am modeling a variable that is not observed but I assume it depends on the observed variable in a linear way and I would like to learn the linear transform from the data.
To clarify things I’ve created a minimal (not) working example. In the code below I would like to learn the value of lag_true based on the observations of value_1 and value_2. To do so I create the matrix A_at and try to use the tensorsolve function to approximate value_1 based on value_2. When I run the code I get the following error:

ValueError: Not enough dimensions on input.
...
Inputs values: [array([[ True]]), array([[-0.91893853]]), array([[-0.5]]), 'not shown', 'not shown', array([[0.5]]), array([[-0.69314718]]), array([[-inf]], dtype=float32)]
Outputs clients: [[Sum{acc_dtype=float64}(sigma > 0)]]


So my questions are:

1. Why is the code below not working ?
2. Is there a better modeling way to achieve what I am trying to do ?

Some comments to give some context on what I am actually trying to do:

1. I assume that the structure of the matrix A_at to depend on the value of lag. Namely, for each column the matrix A_at need to have the lag values below the diagonal equal to 1/lag and the other values equal to 0.
2. In the example I am working on I do not compute value_2 directly from value_1 but rather (through more pymc code) I get an approximation of it. What I would like to, is to learn both the map A_at and the approximation in such a way that the equation A \text{ value_1} = \text{value_2} is ‘‘approximately’’ satisfied.
3. I tried to invert the matrix and compute the matrix-vector multiplication ( as commented in the code ) but this is also not working.

Code:

from aesara.ifelse import ifelse
from aesara.tensor.nlinalg import tensorsolve, tensorinv
import aesara.tensor as at
import numpy as np
import pymc as pm

# create mock dataset
lag_true = 2
n_obs = 50

A_np = np.zeros( ( n_obs, n_obs ) )
for i in range( n_obs ):
A_np[ i:i+lag_true, i ] = 1/lag_true

value_1 = np.random.normal( loc = 0, scale = 1, size = n_obs )
noise = np.random.rand( n_obs ) / 10
value_2 = np.dot( A_np, value_1 ) + noise

# try to learn lag_true with pymc given sell_out and sell_in
with pm.Model() as model:

# random variable for lag
lag = pm.Deterministic( 'lag', 1 + pm.Binomial( 'lag-1', n = 6, p = 0.5 ) )

# data
value_1 = pm.MutableData( 'value_1', value_1 )
value_2 = pm.MutableData( 'sell_2', value_2 )

#create aesara matrix
final_list = []
for i in range( n_obs ):
temp_list = []
for j in range( n_obs ):
if i < j:
temp_list.append( 0 )
else:
temp_list.append( ifelse( at.lt(i, j + lag ), 1/lag, np.float64(0) ) )

final_list.append( temp_list )
A_at = at.stacklists( final_list )
#A_inv = tensorinv( A_at )

value_1_approx = tensorsolve(A_at, value_2)
#value_1_approx = at.dot( A_inv, value_2 )
#value_1_approx = at.tensordot( A_inv, value_2 )

# sample
sigma = pm.Exponential( 'sigma', 2 )
pm.Normal( "obs", mu = value_1_approx, sigma = sigma, observed = value_1 )
trace = pm.sample( 2000, tune = 2000, target_accept = 0.85 )


1 Like