Model Checks, Model Selection, Multivariate Distributions

Exercise 1: Negative Binomial Distribution Derivation¶

In this exercise you will reproduce the derivation of the negative binomial distribution.

a)

The negative binomial distribution gives the probability distribution

over how many $k$ failures are needed before $r$ successes are reached in a binomial experiment, with expectation and variance

Show that by setting $E[k|r,p]\equiv\lambda$ the probability distribution becomes

with $Var[k|r,\lambda]=\lambda(1+\frac{\lambda}{r})$.

Hint: Use the definition of the binomial coefficient via factorials.

No solution, you should check by yourself whether you arrive at the correct result (see lecture slides).

b)

Show that

Hint: See lecture slides.

The Poisson distribution is

For converging to the Poisson distribution for $r\to\infty$, the following properties need to be shown:

• $\displaystyle{\lim_{r\to\infty} \left(1+\frac{\lambda}{r}\right)^r = e^\lambda}$ (known property of the exponential function, see e.g. Wikipedia)

• $\displaystyle{\frac{(r+k-1)!}{(r-1)!} = \frac{1 \cdot 2 \cdot \dots \cdot (r-1) \cdot r \cdot (r+1) \cdot \dots \cdot (r+k-1)}{1 \cdot 2 \cdot \dots \cdot (r-1)} = (r+1) \cdot \dots \cdot (r+k-1)}$

• $\displaystyle{\lim_{r\to\infty} \frac{r \cdot (r+1) \cdot \dots \cdot (r+k-1)}{(r+\lambda)^k)} = \frac{r^k}{r^k} = 1}$

c) Is it a conjugate family?

The Poisson likelihood forms a conjugated family with a Gamma prior distribution. Does the negative binomial distribution also form a conjugated family with a Gamma-distributed prior?

Negative binomial likelihood (use $\alpha$ instead of $r$ because $r$ is already used in prior):

Gamma prior:

Resulting posterior:

This is not equal to the kernel of the Gamma distribution and consequently the negative binomial likelihood and the Gamma prior are not conjugated. We necessarily need to run MCMC to get samples from the posterior.

Exercise 2: More Toilet Paper Data¶

After your preliminary A/B study (previous exercise sheet), you collected a whole year of weekly data (toilet paper weight in kg) in order to give more precise predictions to the facility management so that they can plan their (reduced) budget better. The data are stored in toilet_paper_data.npy (you may load them with np.load()).

a) Simulate the posterior distribution

Like last time, you use a Normal likelihood $y\sim N(\mu,\sigma^{2})$ and assume weak empirical Bayes priors for $\mu$ and $\sigma$: $\mu\sim N(\overline{y},\hat{\sigma})$ and $\sigma\sim Exp(1/\hat{\sigma}),$ where $\overline{y}$ is the empirical mean and $\hat{\sigma}$ the empirical standard deviation of your data $y$.

Simulate the posterior distribution $P(\mu,\sigma|y)$ with PyMC and visualize the marginal distributions $P(\mu|y)$ and $P(\sigma|y)$. What are the 90% highest density intervals for $\mu$ and $\sigma$?

Empirical mean and standard deviation:

Simulation:

Marginal distributions and HDIs:

With 90% plausibility, $\mu$ is between 143-150 and $\sigma$ is between 13 and 18 (values might change slightly after running the simulation with a different seed). After a full year of data, there is quite some (aleatoric) uncertainty remaining!

b) Visualize the Joint Distribution

Use pm.plot_pair() to also visualize the full joint distribution in a contour plot. Does the posterior look more or less (multivariate) normally distributed? Use np.cov() to compute the covariance matrix of the multivariate normal distribution. Are there strong correlations between $\mu$ and $\sigma$ in the posterior? Hint: np.cov() requires a numpy array with each variable (samples for $\mu$ and samples for $\sigma$) as a column.

Looks normally distributed! Note that variances are not equal and you should actually see an ellipse if you turn on equal aspect ratio.

Covariance matrix is almost diagonal! That means there are almost no correlations between $\mu$ and $\sigma$ also in the posterior!

Correlation matrix:

c) Compute samples from the predictive distribution

Compute samples from the predictive distribution using PyMC and use them to compute the root mean squared error (RMSE) and mean absolute error (MAE) averaged over your data and your predictive distribution as done in the lecture.

Reshape predictions for broadcasting with Numpy:

Compute differences:

RMSE (in kg):

MAE (in kg):

RMSE and MAE are quite different, this could indicate the presence of outliers in the data..

d) Perform a Posterior Predictive Check

You are not sure, whether using the normal distribution as likelihood (or sampling probability) was appropriate. To this end, perform a posterior predictive check as done in the lecture. Are the observed data within the range of typical predictions by the model? If not, what tendencies do you observe? What might be the reason for the discrepancy?

The predictions (red) closely follow a normal distribution. However, the center of the observed data distribution is more to the left! There is an outlier present that can be quickly found (using e.g. a strip plot from seaborn):

This outlier appears to pull the estimate for $\mu$ a bit towards the right. This has not much influence on the estimated mean, but quite some on the estimated standard deviation:

Exercise 3: Robust Estimation with Student’s T Distribution¶

(This exercise is a continuation of Exercise 1)

After asking back, the facility management tells you that there was a large conference with many external people in the building on one day (where you have the outlier in your data). You consider two different alternatives:

1. Remove the outlier and re-use the same normal model.

2. Keep the outlier and use a model that is capable of taking outliers into account.

Bayesians often prefer to be subjective in their choice of model but rather like to stay objective in their data. Modifications of data are always dangerous, since they are more difficult to document or reproduce than changes to the model.

a) Visualize the Student’s t distribution

You choose a Student’s t distribution as likelihood, since it is known to be more robust against outliers than the normal distribution. In addition to mean $\mu$ and standard deviation $\sigma$, the Student’s t distribution also has a parameter $\nu$. In statistical t-tests, $\nu$ has typically a non-trivial interpretation as number of degrees of freedom. Here we take a simpler path and use $\nu$ to set the Student’s t distribution in relation with the normal distribution.

Use PreliZ visualisation functionality to compare different Student’s t distributions with the normal distribution. In particular, set $\mu=0$ and $\sigma=1$ and use $\nu\in\{2,3,5,10,20\}$. For what limit of $\nu$ is the Student’s t distribution equal to the normal distribution? In what differs the Student’s t distribution from the normal distribution? Observe in particular also the tails of the different distributions.

Hint: Use the functions pz.StudentT(mu, sigma, nu).plot_pdf() and pz.Normal(mu, sigma).plot_pdf() and make sure the figure has an appropriate size (e.g. figsize=(8,6)) and region of interest (e.g. plt.xlim(-4,4)) so that the differences are clearly visible.

The Student’s $t$ distribution converges towards the normal distribution for large values of $\nu$ (so in theory for $\nu \to \infty$). For smaller values of $\nu$, the Student’s $t$ distribution has fatter tails than the normal distribution.

b) Modeling with Student’s t distribution

Repeat your modelling efforts from Exercise 1 using a Student’s t distribution as likelihood instead of a normal distribution. To this end, use a $Gamma(s=10, r=1.1)$ prior distribution for $\nu$ (that you should visualize beforehand). What are the 90% highest density intervals for $\mu$, $\sigma$ and $\nu$? From looking at $\nu$, do you expect the resulting Student’s t likelihood to have much fatter tails than the previously used normal likelihood?

Prior for $\nu$:

Simulation:

There is slightly less uncertainty in the mean $\mu$ and considerable less uncertainty in $\sigma$ (that was by definition heavily influenced by the outlier)! $\nu$ is with 90% plausibility between 4 and 11, a considerable difference to a normal distribution with fatter tails:

The fatter tails allow to accomodate the outlier there, while the normal distribution tries to bring the outlier more within its center.

c) Posterior predictive checks

Repeat the posterior predictive test that you did in Exercise 1 for your new model and compare it to the old one with normal likelihood. How do they differ?

Student’s $t$ likelihood seems to make predictions quite a bit closer to the data distribution!

d) Compute RMSE and MAE

Compute again RMSE and MAE. What are the differences to Exercise 1 where a normally distributed likelihood was assumed? Were you expecting smaller or larger values compared to Exercise 1? What could be the reason that RMSE and MAE are not as much smaller as you wish them to be?

RMSE:

MAE:

RMSE and MAE are now of almost equal size. However they are only a bit smaller. Choosing a more complex distribution (3 instead of 2 parameters) comes with more uncertainty in the parameters and consequently more predictive uncertainty for the same amount of data. Even though we have reduced the bias of the model (‘wrong’ normal likelihood assumption), we have now increased its variance by introducing an additional parameter.

e) Compare models

Compare both models (normal vs Student’s t likelihood) in terms of expected log predictive density (ELPD). Can you clearly favour one model of the other?

(this warning is justified!!)

In terms of ELPD, the model with the Student’s $t$ likelihood is not significantly better (by the one-standard-error rule) than the model with the normal likelihood. This would probably change with more data and is due to the additional uncertainties (e.g. in the estimation of $\nu$).

Exercise 4: Railway Switch Maintenance¶

For resource planning (personnel and pre-stored material), SBB wants to assess how often their railway switches turn up defective and have to be maintained. To this end, they collected for 10 years the weekly counts of defective railway switches (fictitious data, stored in railway_switch_data.npy).

In this exercise you want to assess whether it is better to use a Poisson or a negative binomial likelihood for your Bayesian model. Since you know that models with more parameters typically pay a price in terms of additional uncertainty, you proceed very carefully.

a) Run posterior simulations

Run posterior simulations for two models:

1. A model with Poisson likelihood and a $Normal(\overline{y}, \hat{\sigma})$ prior on $\lambda$, where $\overline{y}$ is the empirical data mean and $\hat{\sigma}$ the empirical standard deviation.

2. A model with negative binomial likelihood with the same prior on $\lambda$ and a $Gamma(2,0.1)$ prior on $r$.

Quickly visualize the priors used before the simulation—do they make sense for you? Visualize the marginal posterior distributions of the involved parameters and assess their uncertainty. Hint: PyMC does not use exactly the same parameter names for the prior distributions as used here and in the lecture. Make sure you find the correct mappings.

Poisson likelihood model:¶

Negative binomial likelihood model:¶

Again the joint posterior distribution looks normally distributed!

b) Perform a posterior predictive check

Perform a posterior predictive check of both models and plot the side-by-side. Which model do you prefer after the posterior predictive check?

Posterior predictions:

Posterior predictive checks:

The negative binomial likelihood seems much better suited to model the data!

c) Compute Bayes factors

Compute the Bayes factor between the marginal likelihoods for the two models. Which model explains the data better in your belief? Also try to compute the marginal likelihoods $p(d|\mathcal{M}_{Pois})$ and $p(d|\mathcal{M}_{NegBin})$. Would you expect that they are so extremely small?

Hint: Use SMC as introduced in the lecture and the companion notebook.

It is beyond any discussion that the negative binomial model is better suited (BF=10⁸²..)

The marginal probabilities are extremely small!

This is the reason why we work with log probabilities, to at least being able to evaluate the Bayes factor. The marginal probabilities are so small because they are the product of the probabilities of all data points under the posterior distribution:

Imagine each data point has on average 80% probability under the model (which is high). Then the resulting marginal probability is

extremely small values are to be expected!

d) Compare models

Compare the models in terms of ELPD. Is there a clear winner when involving the one-standard-error rule?

Clearly it is recommended to use the negative binomial model! This time we have enough data to reduce the uncertainty of the additional parameter.

e) Compute predictive RMSE and MAE

Compute predictive RMSE and MAE for both models and compare them.

RMSE:

MAE:

Both negative binomial RMSE and MAE are a bit smaller!

Exercise 5: Estimating Churn Proportions¶

You work at a national health insurance company. Recently, significantly more customers than usual have terminated their insurance policies with your company and changed to somebody else. In the insurance jargon, the number of customers terminating their contract is called churn. The management has ordered an investigation that was carried out by the sales department: they have tried to contact 200 customers that have left and asked them about the reason why they had terminated their contract. They could actually reach 150 customers and out of these 112 gave a suitable answer. Finally, the sales departement has put the different churn reasons into three main categories:

1. premium (‘prämie’) too high ($k_{1}=82$ cases)

2. unsufficient coverage ($k_{2}=24$ cases)

3. unsuitable insurance models ($k_{3}=6$ cases)

As a data scientist, you are aware that this is a restricted sample and decide to use a Bayesian model with multinomial likelihood and flat Dirichlet prior to model the true underlying proportions $\pi_{1}$, $\pi_{2}$ and $\pi_{3}$ of your entire insurant population including uncertainty.

a) Simulate the posterior

Simulate the posterior with PyMC and visualize the marginal posterior distributions of $\pi_{1}$, $\pi_{2}$ and $\pi_{3}$. What 90% HDIs can you communicate back to the management? Formulate a sentence.

I believe to 90% that for around 67-78% of the churn customers the premium is too high, for around 17-27% the coverage was unsufficient and around 0.03-0.09% found our insurance models unsuitable.

b) Use the Dirichlet-multinomial conjugacy

Use Dirichlet-multinomial conjugacy and Preliz to visualize your posterior on a probability simplex.

c) Switch to a different prior

Your superior is unsure about your results. She says: “It’s very typical that people switch contracts because of the premium. I would have used a Dirichlet ([10,1,1]) prior that reflects this.”

Use Preliz to visualize this prior on a simplex and marginalized (to look at $\pi_{1}$ specifically) and compute the prior mean and standard deviation for $\pi_{1}$. Finally, re-run your simulation and show that your results are still pretty similar and that your posterior is quite insensitive to this choice of prior.

90% HDIs are only changed by a few percent! This will not have a big impact on my communication.

Exercise 6: The Multivariate Normal Distribution¶

Even though the multivariate normal distribution plays a big role in Bayesian inference, it will not play a major role in the first part of this lecture, since the bambi library does everything for us in hiding complexity. However, the multivariate normal distribution will become more important, when you will work with Gaussian processes in the second half.

The following three exercises are meant to get to know the multivariate normal distribution a little bit better.

a) The multivariate normal distribution

Load the credit rating data (from credit_data.csv in the exercise materials) and select ‘Income’ as $x$ and ‘Rating’ as $y$. Compute (classical) estimators for $\mu=(\begin{matrix}\mu_{x}\\ \mu_{y}\end{matrix})$ and the covariance matrix $\Sigma$. Visualize the resulting normal distribution $N(\hat{\mu},\hat{\Sigma})$ using a contour plot and overplot with the data. Ponder (qualitatively) whether the normal distribution is a good model for this data.

Hints: You may use np.cov() to compute the covariance and code in the notebook ‘week5_multivariate_distributions.ipynb’ to produce the contour plot.

Compute means and covariance:

Plot contour and data (removed equal axis dimensions to make the plot a bit more appealing):

b)

Show that if $x_{1}\sim N(\mu_{1},\sigma_{1}^{2})$ and $x_{2}\sim N(\mu_{2},\sigma_{2}^{2})$, then

Hints: Compute $E[x_{1}+x_{2}]$ and $Var[x_{1}+x_{2}]$ and use linearity in the first and that $Var[x_{1}+x_{2}]=Var[x_{1}]+Var[x_{2}]$ for independent $x_{1}$ and $x_{2}$.

c)

Let $X\sim N(\mu,\Sigma)$ be a random vector. Show that $AX+b\sim N(A\mu+b,A\Sigma A^{T}).$

Hints: Compute $E[A\mu+b]$ and use linearity. Compute $Cov(A\mu+b)$ and use that $Cov(A\mu+b)=Cov(A\mu)$ and $\Sigma=Cov(X)=E[(X-E[X])(X-E[X])^{T}]=E[XX^{T}]-E[X]E[X]^{T}$.