Hi Jordan
A few small issues that I can see. The shape of y should be 1dimensional so (150,) likewise beta should have shape=(X.shape[1]) rather than shape=(X.shape[1], 1)
A bigger issue is the width of the priors you use with Normal(..., sd=10). I would recommend sd=1. Given your inputs, X, a beta value of more than 10 doesn’t make sense, apriori we’d think a value in range abs 0-5 might make sense. So a Normal(..., sd=1) can cover that whereas sd=10 is too wide.
.
The sampling should be much better after these changes and you can lower the number of tuning samples to a more reasonable number like 1,000, rather than 10,000