Question on how to model spike and slab priors

Yes, that is what I meant.
I see youve updated your post with your model.

What are your sampling options? Because you cannot use NUTS for discrete variables, you need to (dramatically) increase your tuning steps. When I choose

trace = pm.sample(1000, tune=12000, nuts={"target_accept": 0.99})

with (this is adapted with your code)

sigma_1 = pm.HalfNormal('sigma_1',10, shape=shape)
pi = pm.Beta('pi',2,7)
gamma = pm.Bernoulli('gamma',p = pi,shape = shape)
mixture_sd = pm.Deterministic('mixture_sd', pm.math.switch(gamma > 0.5, sigma_1, 0.001))
beta = pm.Normal('beta',mu = 0,sigma = mixture_sd, shape = shape)

I do get the right beta coefficients with my toy dataset.

And I think the next step would be to not use match.switch, but replace it with a continous function like tt.net.sigmoid, to improve sampling, see Frequently Asked Questions