Sure, it’s like this:
with pm.Model() as model1:
bounded_laplace = pm.Bound(pm.Laplace, lower=0.0)
alpha = bounded_laplace('alpha', mu = 51704326, b = 767593)
beta_param = SSBeta(name = 'beta_testing',
alpha = 1.578,
beta = 0.855,
loc = 0.0004,
scale = 7.338e-5)
pred = pm.Gamma("test_var", alpha = alpha, beta = beta_param,
observed=np.array([1.72e6, 1.58e6, 1.64e6, 1.59e6]))
trace = pm.sample(400)