Shape error when making out-of-sample predictions

Questions version agnostic
,
7

I tried in a simplified scenario:

import pymc as pm
import numpy as np

N = 2
N_test = 3
D = 4

X = np.random.normal(0,1, size=(N,D))
y = np.random.normal(0,1, size=(N,))
X_test = np.random.normal(0,1, size=(N_test,D))

model = pm.Model()
with model:
    X = pm.MutableData ('X', X)
    N = pm.MutableData('N',N)
    N_int = np.int64(N.get_value())
    factors = pm.Normal('factor',mu=0, sigma = 1, shape = (N_int,))
    output = np.ones((N_int,))*factors
    output = pm.Deterministic('output',output)
    likelihood = pm.Normal('y', mu=output, sigma=1, observed=y, shape=output.shape)
    
    if __name__ == "__main__":
        trace = pm.sample()
        pm.set_data({'X': X_test, 'N':np.shape(X_test)[0]})
        samples = pm.sample_posterior_predictive(trace,predictions=True)
        print(samples.predictions)

Print of samples.predictions looks like

Dimensions:  (chain: 4, draw: 1000, y_dim_2: 2)
Coordinates:
  * chain    (chain) int32 0 1 2 3
  * draw     (draw) int32 0 1 2 3 4 5 6 7 8 ... 992 993 994 995 996 997 998 999
  * y_dim_2  (y_dim_2) int32 0 1
Data variables:
    y        (chain, draw, y_dim_2) float64 1.394 -2.583 ... -1.247 -0.0293

which is the wrong shape (it should be 3 not 2). The problem seems to be that I can’t make predictions with a model which is an explicit function of the shape of the data.