Q52. Bayesian logistic regression

Topic: Markov chains, Gibbs sampling, and Metropolis–Hastings · more questions on this topic
‹ Q51 Q59 ›

[12 marks]

A researcher models whether a manufactured component fails a stress test. For component i, let y_i=1 indicate failure and let x_i denote its centred test temperature. Conditional on \boldsymbol\beta=(\beta_0,\beta_1)^\top, assume Y_i\mid\boldsymbol\beta\overset{\mathrm{ind}}{\sim} \operatorname{Bernoulli}(p_i),\qquad p_i=\frac1{1+\exp[-(\beta_0+\beta_1x_i)]},\quad i=1,\ldots,n. Assign independent priors \beta_0\sim N(0,4),\qquad\beta_1\sim N(0,9), where the second parameter denotes the variance. Treat the x_i as fixed.

(a) Write down the likelihood L(\boldsymbol\beta;\mathbf y,{\boldsymbol{x}}). [2 marks]

(b) Derive the posterior density of \boldsymbol\beta, up to proportionality. Hence obtain its log density up to an additive constant, expressed explicitly in terms of \beta_0,\beta_1,x_i,y_i. [3 marks]

A random-walk Metropolis–Hastings sampler proposes \boldsymbol\beta^*\mid\boldsymbol\beta^{(t)} \sim N_2\!\left(\boldsymbol\beta^{(t)},s^2I_2\right), where s>0 is fixed and I_2 is the 2\times2 identity matrix.

(a) Starting from the general Metropolis–Hastings acceptance probability, simplify it for this proposal. Explain any cancellations and express your answer using the log posterior from part (b). [3 marks]

(b) Give pseudocode for the sampler. Specify what happens after both an accepted proposal and a rejected proposal. [2 marks]

(c) Explain how choosing s extremely small or extremely large can impair the sampler’s efficiency. [2 marks]

No implementation or numerical simulation is required.

Worked answer (click to show)

Definitions. Write \boldsymbol\beta=(\beta_0,\beta_1)^\top, \eta_i=\beta_0+\beta_1x_i, and p_i=e^{\eta_i}/(1+e^{\eta_i}). A Bernoulli observation contributes p_i^{y_i}(1-p_i)^{1-y_i} to the likelihood. The fixed x_i are covariates, the binary outcomes are conditionally independent, and the independent Normal priors have variances 4 and 9. I_2 is the identity matrix; the proposal covariance is fixed at s^2I_2.

The logistic probabilities satisfy p_i=\frac{e^{\eta_i}}{1+e^{\eta_i}},\qquad 1-p_i=\frac1{1+e^{\eta_i}}.

(a) Likelihood. For observed y_i\in\{0,1\}, conditional independence gives \boxed{ L(\boldsymbol\beta;\mathbf y,{\boldsymbol{x}}) =\prod_{i=1}^n p_i^{y_i}(1-p_i)^{1-y_i} =\prod_{i=1}^n\frac{\exp\{y_i(\beta_0+\beta_1x_i)\}} {1+\exp(\beta_0+\beta_1x_i)}. } The covariates are fixed, so no probability density for the x_i is included in this likelihood.

(b) Posterior and log posterior. Independence of the priors gives p(\boldsymbol\beta)=p(\beta_0)p(\beta_1) \propto\exp\left(-\frac{\beta_0^2}{2(4)} -\frac{\beta_1^2}{2(9)}\right) =\exp\left(-\frac{\beta_0^2}{8}-\frac{\beta_1^2}{18}\right). The denominators use the stated variances 4 and 9, not their squares. Multiplying this prior kernel by the likelihood yields \boxed{ \pi(\boldsymbol\beta\mid\mathbf y,{\boldsymbol{x}}) \propto \left[\prod_{i=1}^n \frac{\exp\{y_i(\beta_0+\beta_1x_i)\}} {1+\exp(\beta_0+\beta_1x_i)}\right] \exp\left(-\frac{\beta_0^2}{8}-\frac{\beta_1^2}{18}\right), } with support \boldsymbol\beta\in\mathbb R^2. Taking logarithms, \log\pi(\boldsymbol\beta\mid\mathbf y,{\boldsymbol{x}}) =h(\boldsymbol\beta)+C, where C does not depend on \boldsymbol\beta and \boxed{ h(\boldsymbol\beta) =\sum_{i=1}^n\left[ y_i(\beta_0+\beta_1x_i) -\log\{1+\exp(\beta_0+\beta_1x_i)\}\right] -\frac{\beta_0^2}{8}-\frac{\beta_1^2}{18}. } Thus h is the log posterior up to an additive constant; it includes the log likelihood and the log prior.

(c) Acceptance probability. For current state \boldsymbol\beta and proposed state \boldsymbol\beta^*, the general rule is A(\boldsymbol\beta,\boldsymbol\beta^*) =\min\left\{1, \frac{\pi(\boldsymbol\beta^*\mid\mathbf y,{\boldsymbol{x}}) q(\boldsymbol\beta\mid\boldsymbol\beta^*)} {\pi(\boldsymbol\beta\mid\mathbf y,{\boldsymbol{x}}) q(\boldsymbol\beta^*\mid\boldsymbol\beta)}\right\}. The bivariate Normal random-walk proposal has density q(\boldsymbol\beta^*\mid\boldsymbol\beta) =\frac1{2\pi s^2} \exp\left[-\frac{\|\boldsymbol\beta^*-\boldsymbol\beta\|^2}{2s^2}\right]. Since the squared Euclidean distance is unchanged when the two states are exchanged, and s^2I_2 is fixed, q(\boldsymbol\beta^*\mid\boldsymbol\beta) =q(\boldsymbol\beta\mid\boldsymbol\beta^*). The proposal-density terms therefore cancel. The posterior’s normalising constant also cancels, giving \boxed{ A(\boldsymbol\beta,\boldsymbol\beta^*) =\min\{1,\exp[h(\boldsymbol\beta^*)-h(\boldsymbol\beta)]\}. }

(d) Pseudocode. Start from any finite \boldsymbol\beta^{(0)}\in\mathbb R^2, for example (0,0)^\top, and choose fixed s>0. For t=0,\ldots,T-1:

1. Draw \boldsymbol\beta^*\sim N_2(\boldsymbol\beta^{(t)},s^2I_2).

2. Compute \Delta=h(\boldsymbol\beta^*)-h(\boldsymbol\beta^{(t)}).

3. Independently draw U\sim\operatorname{Uniform}(0,1).

4. If \log U\le\min\{0,\Delta\}, accept and set \boldsymbol\beta^{(t+1)}=\boldsymbol\beta^*. Otherwise reject and set \boldsymbol\beta^{(t+1)}=\boldsymbol\beta^{(t)}.

5. Store \boldsymbol\beta^{(t+1)}, including repeated states following rejection.

The log-scale comparison avoids exponentiating \Delta. Keeping only accepted proposals would not give the intended sequence of posterior Markov-chain states. Allow for convergence from the initial state before using the output for posterior summaries.

(e) Effect of proposal scale.

• Extremely small s: proposals tend to be accepted often, but move only a short distance. The chain explores the posterior slowly and successive states are highly correlated, so many iterations provide little additional information.

• Extremely large s: proposals often land in regions of much lower posterior density and are rejected. The chain then repeats its current state many times, again exploring slowly.

Both extremes reduce efficiency. A high acceptance rate alone is not evidence of efficient sampling; useful movement and mixing also matter. For any fixed s>0, the algorithm still targets the same posterior—the choice changes how efficiently it explores that distribution.