Model with a loop - scan conversion

Questions v5
2

Hey hey

You’re on the right track, but your mental model of what pytensor is is wrong. Amusingly, my very first post on this forum made the same mistake, so you’re in good company.

Pytensor works by making a graph of symbolic computation. When you do e.g. z = x + y, z is not a number, it’s a set of instructions describing how x, y, and z are related. You build a computation up bit by bit, using and re-using these computational graphs. When you call pytensor.function, you signal that you’re done build and you want to do actual numeric computation from the inputs to the outputs.

So, your proximate error is caused because you are passing symbolic values like alpha to your compiled function, which is done with symbolic stuff and is now ready for numeric values.

But the ultimately error is that you shouldn’t be doing anything with pytensor.function in the context of a pymc model. PyMC is going to compile some function itself, and it can’t see into or use your compiled function. But you don’t need this – instead, you need to keep everything symbolic from start to finish. For example:

import pytensor 
import pytensor.tensor as pt
import numpy as np
import pymc as pm

n = 100
temp_data = np.random.normal(size=(n, ))

with pm.Model() as model:
    alpha = pm.Normal("alpha", mu=2.5, sigma=0.5)
    t_eq  = pm.Normal("t_eq", mu=-1.5, sigma=0.5)
    s0    = pm.Normal("s0", mu=100, sigma=20)
    rho0  = pm.TruncatedNormal("rho0", mu=0.5, sigma=0.3, lower=-1, upper=1)
    sigma = pm.TruncatedNormal("sigma", mu=8, sigma=3, lower=0, upper=15)

    temps = pm.Data('temps', temp_data)

    def scan_step(previous_level, alpha, temps, t_eq):
        return previous_level + alpha * (temps - t_eq)
    
    sea_levels = pytensor.scan(scan_step,
                               outputs_info = s0,
                               sequences = temps,
                               non_sequences = [alpha, t_eq],
                               return_updates=False)
    
    hay_stds_diag = pt.eye(n)
    
    sigma_ar_1 = sigma **2 / (1 - rho0 ** 2)
    powers = pt.arange(1, n+1)
    cov_matrix_hay = sigma_ar_1 * rho0 ** pt.abs(pt.sub.outer(powers, powers)) + hay_stds_diag
    
    sea_level_hat = pm.MvNormal(name='sea_levels',
                                mu=sea_levels,
                                cov=cov_matrix_hay)

Some other issues though:

You seem to be computing an AR(1) model twice, once in the mean, then again in the covariance. It seems like you’re mixing a model defined once as a generative model (the AR(1) scan) and a likelihood (the MvNormal with AR(1) covariance matrix)

Using a big covariance matrix to define the observed variable will end up being really slow compared to using pm.AR. You’re essentially making a Gaussian Process by hand. Using a GP would give you access to approximations that could speed things up (like HSGP)

Nitpick: What you call “likelihood” isn’t a likelihood. It’s a symbolic random variable that generates draws from your model. You can check that pm.draw(likelihood) does not give back evaluations of the likelihood function. This distribution will be used to compute the likelihood term during estimation, but it’s not there yet.