Adding constraint to the parameter

You can use potential sure but sampling will be really inefficient.

pm.Potential(‘p_min_potential’, tt.switch(tt.eq(tt.sum(p), 1.), 0, -np.inf))

You can also normalized it p = p/p.sum() but Dirichlet is the best solution.