Send your solutions to: vavraj@karlin.mff.cuni.cz by 3 November, 9:00
When measuring the speed of \(n = 50\) cars on Rohanské nábřeží, opposite of the MFF building, we obtained the sample mean \(\overline y = 59.6\) and the sample standard deviation of \(s = 11.1\) km/h (the measurements were taken in the period before installation of two traffic lights at the corners of the River Garden office building).
Examine the posterior distribution of the mean (expected value) of the speed of cars at Rohanské nábřeží while assuming a normal distribution \(\mathsf{N}(\mu,\sigma^2)\) of the speed. As a prior for \(\mu\) and \(\tau = \sigma^{-2}\) consider a (semi-conjugate) distribution: \[ p(\mu,\,\tau) = p(\mu)\,p(\tau), \] where \(p(\mu) \propto 1\) and \(\tau\sim\mathsf{\Gamma}(g,\,h)\), \(p(\tau) \propto \tau^{g-1}\exp(-\,h\,\tau)\). For \(g\) and \(h\) consider the following choices:
As a check, calculate the ordinary (i.e., non-bayesian, frequentist) 95% confidence interval for the unknown mean \(\mu\) and also for the standard deviation \(\sigma\) and inverse variance \(\tau\).
We can derive the CI from the following distributions: \[ \overline{Y}_n \sim \mathsf{N}\left(\mu,\, \frac{\sigma^2}{n}\right) \qquad \text{and} \qquad \frac{(n-1)S_n^2}{\sigma^2} \sim \chi_{n-1}^2. \]
Derive the marginal posterior distribution for \(\mu\) and plot its density (in one figure) for the above choices \(g\) and \(h\). What is the posterior probability that the mean of the speed is greater than 55 km/h (again for three different choices of \(g\) and \(h\))?
We need to marginalize the following pdf with respect to \(\tau\):
\[ p(\mu, \tau | \mathbb{Y}) \propto p(\mathbb{Y} | \mu, \tau) p(\mu, \tau) \propto \tau^{\frac{n}{2}} \exp \left\{-\frac{\tau}{2}\sum\limits_{i=1}^n(Y_i - \mu)^2\right\} \cdot 1 \cdot \tau^{g-1} \exp \left\{-\tau h\right\}. \] This envolves integration \(\int\limits_{0}^\infty \cdots \mathrm{d}\tau\), which can be solved by using knowledge of \(\Gamma(\alpha, \lambda)\) distribution.
Hint: Up to a multiplicative constant, you should recognize the remaining function of \(\mu\) as pdf of Location-scale Student distribution.
Derive the marginal posterior distribution for \(\tau\) and plot its density (in one figure) for the above choices \(g\) and \(h\). In the second plot, draw the marginal posterior densities for \(\sigma = \sqrt{1/\tau}\). Note that it is not necessary to explicitly derive the formula for the posterior density of \(\sigma\) for the purpose of plotting once you have a computer function to calculate the functional values of the posterior density of \(\tau\).
We need to marginalize the following pdf with respect to \(\mu\):
\[ p(\mu, \tau | \mathbb{Y}) \propto p(\mathbb{Y} | \mu, \tau) p(\mu, \tau) \propto \tau^{\frac{n}{2}+g-1} \exp \left\{-\tau \left[\frac{1}{2}\sum\limits_{i=1}^n(Y_i - \mu)^2 + h\right]\right\}. \] This envolves integration \(\int\limits_{-\infty}^\infty \cdots \mathrm{d}\mu\), which can be solved by using knowledge of \(\mathsf{N}(\mu,\sigma^2)\) distribution.
Hints: Up to a multiplicative constant, you should recognize the remaining function of \(\tau\) as pdf of some well known distribution. To obtain the marginal posterior density of \(\sigma\) use transformation theorem.
Draw (in three different plots) plot of a joint posterior density of \((\mu,\,\tau)\) for three different choices of \(g\) and \(h\). Note that you can use the relation \(p(\mu,\,\tau\,|\,\mathbb{Y}) = p(\mu,\,\tau\,|\,\mathbb{Y})\,p(\tau\,|\,\mathbb{Y})\) to calculate the functional values of \(p(\mu,\,\,\tau\,|\,\mathbb{Y})\).
Hint: Use the following functions:
x; y # grids on x and y axis
z <- outer(x, y, fun) # fun = function of x and y to be evaluated
image(x, y, z) # creates colourful plot
contour(x, y, z, add = T) # adds contours
# alternatively
persp(x, y, z) # 3D plotDraw (in three different plots) plot of a joint posterior density of \((\mu,\,\sigma)\) for three different choices of \(g\) and \(h\). Note that you can use the relation \(p(\mu,\,\sigma\,|\,\mathbb{Y}) = p(\mu,\,\sigma\,|\,\mathbb{Y})\,p(\sigma\,|\,\mathbb{Y})\) to calculate the functional values of \(p(\mu,\,\,\sigma\,|\,\mathbb{Y})\). Moreover, \(p(\mu\,|\,\sigma,\,\mathbb{Y}) = p\bigl(\mu\,\big|\,\sigma^{-2},\,\mathbb{Y}\bigr)\).
Calculate the 95% ET (equal-tail) credible interval for \(\mu\), \(\tau\), \(\sigma\) and each choice of \(g\) and \(h\).
Calculate (numerically, if needed) the 95% HPD (highest posterior density) credible interval for \(\mu\), \(\tau\), \(\sigma\) and each choice of \(g\) and \(h\). Do those intervals differ from the ET credible intervals?
Write a short function which can be used to generate pseudorandom numbers from the joint posterior distribution \((\mu,\,\tau)\) (specify the prior hyperparameters as arguments of the function). Note once again that \(p(\mu,\,\tau\,|\,\mathbb{Y}) = p(\mu\,|\,\tau,\,\mathbb{Y})\,p(\tau\,|\,\mathbb{Y})\).
Based on the simulated sample from the posterior distribution (of a length at least 10,000), calculate again (now Monte Carlo estimates of) the 95% ET credible intervals of parameters \(\mu\), \(\tau\) and \(\sigma\) (again for the three choices of hyperparameters \(g\) and \(h\)). Do those intervals differ from the intervals calculated in point 6?
Extend the above function to generate also from the posterior predictive distribution of the speed \(Y_{n+1}\) of another car: \[\begin{equation}\label{eqpredhEN} p(y_{n+1}\,|\,\mathbb{Y}) = \int p(y_{n+1},\,\mu,\,\tau\,|\,\mathbb{Y})\,\mathrm{d}(\mu,\,\tau) = \int p(y_{n+1}\,|\,\mu,\,\tau,\,\mathbb{Y})\,p(\mu,\,\tau\,|\,\mathbb{Y})\,\mathrm{d}(\mu,\,\tau) = \int p(y_{n+1}\,|\,\mu,\,\tau)\,p(\mu,\,\tau\,|\,\mathbb{Y})\,\mathrm{d}(\mu,\,\tau). \end{equation}\] Consequently, calculate (based on at least 10,000 simulated values from the predictive distribution \(p(y_{n+1}\,|\,\mathbb{Y})\)) the 95% ET credible intervals for \(Y_{n+1}\) (again, while considering the three choices of hyperparameters \(g\) and \(h\)). How would you interpret the credible intervals for \(Y_{n+1}\)?
Remind that to generate from the distribution \(p(y_{n+1}\,|\,\mathbb{Y})\) it might be useful to generate from seemingly a more complicated joint distribution \[\begin{equation*} p(y_{n+1},\,\mu,\,\tau\,|\,\mathbb{Y}) = p(y_{n+1}\,|\,\mu,\,\tau,\,\mathbb{Y})\,p(\mu,\,\tau\,|\,\mathbb{Y}) = p(y_{n+1}\,|\,\mu,\,\tau)\,p(\mu,\,\tau\,|\,\mathbb{Y}). \end{equation*}\]
Use the simulation and approximate values of the predictive density \(p(y_{n+1}\,|\,\mathbb{Y})\) (again for the three choices of hyperparameters \(g\) and \(h\)) and draw them in one plot. Note that the MC estimate of the predictive density can be obtained while using the relationship from Task 9.
Extend a bit more the function which generates from the posterior distribution \(p(\mu,\,\tau\,|\,\mathbb{Y})\) and calculate the Monte Carlo estimate of a probability that another car exceeds the allowed speed by more than 30 km/h (and the driver loses their driving license for certain period of time).
Note an analogue of the relationship from Task 9: \[\begin{equation*} \mathsf{P}(Y_{n+1} > 80\,|\,\mathbb{Y}) = \int \mathsf{P}(Y_{n+1} > 80\,|\,\mu,\,\tau,\,\mathbb{Y})\, p(\mu,\,\tau\,|\,\mathbb{Y})\, \mathrm{d}(\mu,\,\tau) = \int \mathsf{P}(Y_{n+1} > 80\,|\,\mu,\,\tau)\,p(\mu,\,\tau\,|\,\mathbb{Y})\,\mathrm{d}(\mu,\,\tau). \end{equation*}\]