How to Get Predictions from Your Fitted Bayesian Model in Python + R

Last week we built our first Bayesian linear regression model using Stan. This week we continue using the same model and data set from the Spotify API to generate and visualise our predictions.

Once we have our model we can generate new predictions. With a Bayesian model we don't just get a prediction but a population of predictions.

Which we can visualise as a distribution:

Which yields the plot you see in the cover image. Now we will replicate this process using PyStan in Python.

You can find the definition of the stan_code and data in last weeks edition of Data Science Code in Python + R. Note that we are taking a different approach to generate our predictions that will require you to import numpy to execute.

Now that we have our predictions we will create a data frame and build our plot using seaborn which you will need to import as well to run the code below.

Which generates the following, somewhat narrower distributions for the popularity. I haven't yet figured out why the distributions look quite different so if you think you know why please let me know in the comments.

In the R code you will note some code that has been commented out. This takes longer to run but sets the priors to be equivalent to the Stan code we are using in the Python script. The end result was not much different from when using the default priors set by the RStanArm package.

Next week we will continue with our dive into running Stan in Python and R with how to construct and visualise credible intervals for our model and one way to select the best fitting model.

Until then, if you are found this post helpful like and subscribe to Data Science Code in Python + R and let me know what you would like to see in future posts in the comments.

~Matt