Cannot convert to variable (Aesara scan)

Questions
1

Below is a MWE – I know I could either vectorize the model or use the AR class directly, but I’m trying to use this as

(1) a chance to learn more about some of the internals/features of aesara
(2) the baseline for a more complicated computation (that won’t necessarily allow me to get away with using built-in features or vectorizing).

import aesara
import aesara.tensor as aet
import numpy as np
import pymc3 as pm

ar_data = np.array(
    [
        0., 0., 0.06535454, 0.04804798, -0.00423464,
        0.00828503, -0.06116375, -0.09768569, 0.00700567
    ]
)

ar_bh = pm.Model()


with ar_bh:
    # Data
    data_bh = pm.Data("data", ar_data)

    # Parameters
    rho_1 = pm.Uniform("rho_1", -1, 1)
    rho_2 = pm.Uniform("rho_2", -1, 1)
    sigma = pm.HalfNormal("sigma", 3)

    # Update via scan
    def one_step_mu(ym1, ym2, rho_1, rho_2):
        return rho_1*ym1 + rho_2*ym2

    mu = aet.dvector("mu")
    _mu = aesara.scan(
        one_step_mu,
        sequences=[dict(input=data_bh, taps=[-1, -2])],
        outputs_info=[mu],
        non_sequences=[rho_1, rho_2]
    )

    obs = pm.Normal("obs", mu, sigma, observed=data_bh)

This code produces a NotImplementedError: Cannot convert rho_1 ~ Uniform to a tensor variable.

Is this expected? Can random variables not be an input to scan?

(Similar code without being in a model context seems to work)

def one_step(ym1, ym2, rho_1, rho_2):
    return rho_1*ym1 + rho_2*ym2

y = aet.vector("y")
rho_1 = aet.dscalar()
rho_2 = aet.dscalar()

results, updates = aesara.scan(
    one_step, sequences=[dict(input=y, taps=[-1, -2])],
    non_sequences=[rho_1, rho_2]
)
compute_mu = aesara.function(inputs=[y, rho_1, rho_2], outputs=results, updates=updates)
mu = compute_mu(ar_data, 0.9, -0.2)
2

With some slight modifications, the code seems to work when called with theano rather than aesara.

import theano as aesara  # Forgive the naming abuse -- It keeps code the same
import theano.tensor as aet  # ^See above...
import numpy as np
import pymc3 as pm

ar_data = np.array(
    [
        0., 0., 0.06535454, 0.04804798, -0.00423464,
        0.00828503, -0.06116375, -0.09768569, 0.00700567
    ]
)
ar_bh = pm.Model()


with ar_bh:
    # Data
    data_bh = pm.Data("data", ar_data)

    # Parameters
    rho_1 = pm.Uniform("rho_1", -1, 1)
    rho_2 = pm.Uniform("rho_2", -1, 1)
    sigma = pm.HalfNormal("sigma", 3)

    # Update via scan
    def one_step_mu(ytm1, ytm2, rho_1, rho_2):
        return rho_1*ytm1 + rho_2*ytm2

    _mu, update = aesara.scan(
        one_step_mu,
        sequences=[dict(input=data_bh, taps=[-1, -2])],
        non_sequences=[rho_1, rho_2]
    )

    obs = pm.Normal("obs", _mu, sigma, observed=data_bh[2:])