Perhaps I’m misunderstanding, but it seems there isn’t any randomness in this setup. That is, given that you know the equation y= X\beta, where y is the channel A spend, X is the channel_spend_A_data, and \beta is channel_A_coef, you can deterministically solve the linear system for \beta = X^{-1}y. In code:
import pytensor.tensor as pt
channel_A_coef = pt.linalg.inv(channel_A_spend) @ channel_A_uplift
That would remove the constraint entirely. You could actually just pre-compute this number using numpy outside of your model to save your CPU from doing the matrix inversion on each logp call, since it’s totally deterministic.
Given the shapes of your dummy data though, I guess you can’t actually invert channel_A_spend (it’s just a column vector?). So then you could obtain channel_A_coef via minimization and use that:
from scipy.optimize import minimize_scalar
def objective(beta):
return (beta * np.array(channel_A_spend) - np.array(channel_A_uplift)).sum() ** 2
res = minimize_scalar(func, bounds=[1e-4, 100])
Which you could then deterministically plug into the model, and estimate everything else.