Simple way to making a prediction using PyMC model

Questions
4

Sure. Or you could just sum across prices, but keep the sum separate for each sample in your posterior:

posterior_rev = np.sum([p * (mu * price) for price in test['Price']], axis=0)

Then I would probably feed this to arviz’s plot_dist().