Learn linear transformation with pymc

Questions v5
1

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