Metropolis-Hastings for linear regression, prior on sigma?

Question

For the sake of curiosity, I'm trying to build a Metropolis-Hastings sampler for the purposes of Bayesian linear regression. Below, you'll note my script and more specifically, in-line comments that noting to comment in/out various lines in order to change the behavior of the script.

As is, the sampler iteratively proposes a change to either b (the slope), or a (the slope intercept.) However, these parameters of the linear function y=bx+a, are not updated simultaneously. It works great! However, as is, the script does not propose changes to s, the standard deviation of the linear function. When I alter the code block as detailed, the sampler fails to change. It simply stagnates at the initial possible value.

My questions are:

(1) What are the benefits of sampling different sigma values? I get a pretty good understanding of b and a with s constant.
(2) Am I proposing changes to s wrong? I understand that it cannot be negative, but it also needs to be sampled from a symmetric distribution, allowing for increases and decreases. I've used absolute value of the current value plus some random change. (-0.15 -> 0.15) (3) Is there a better prior choice for sigma? I'm using inverse gamma. Also, you'll note that none of my distribution functions involve normalizing constants as this generally isn't necessary in MH.

My code:

import numpy as np
import random

def normalPDF(x,mu,sigma):
  num = np.exp((x-mu)**2/-2*sigma**2)
  return num

def invGamma(x,a,b):
  non_zero = int(x>=0)
  func = x**(a-1)*np.exp(-x/b)
  return non_zero*func

def lr_mcmc(X,Y,hops=10_000):
  samples = []
  
  curr_b = 1
  curr_a = 1
  curr_s = 1

  prior_b_curr = normalPDF(x=curr_b,mu=2,sigma=1)
  prior_a_curr = normalPDF(x=curr_a,mu=1,sigma=1)
  prior_s_curr = invGamma(x=curr_s, a=2,b=2)
  
  log_lik_curr = sum([np.log(normalPDF(x=curr_b*x + curr_a,mu=y,sigma=curr_s)) for x,y in zip(X,Y)])
  current_numerator =  log_lik_curr + np.log(prior_a_curr) + np.log(prior_b_curr) + np.log(prior_s_curr)

  count = 0
  for i in range(hops):
    samples.append((curr_b,curr_a,curr_s))    

    if count == 0:
      mov_b = curr_b + random.uniform(-0.25,0.25)
      mov_a = curr_a 
      mov_s = curr_s
      count += 1

    elif count == 1:
      mov_a = curr_a + random.uniform(-0.25,0.25)
      mov_b = curr_b
      mov_s = curr_s

      # to change behavior:
      # count += 1 # uncomment line 
      count = 0    # comment line out

      # to change behavior, uncomment below code block:
    # else:
    #   mov_s = np.abs(curr_s + random.uniform(-0.25,0.25))
    #   mov_b = curr_b
    #   mov_a = curr_a
    #   count = 0

    prior_b_mov = normalPDF(x=mov_b,mu=2,sigma=1)
    prior_a_mov = normalPDF(x=mov_a,mu=1,sigma=1)
    prior_s_mov = invGamma(x=mov_s,a=2,b=2)
    log_lik_mov = sum([np.log(normalPDF(x=mov_b*x + mov_a,mu=y,sigma=mov_s)) for x,y in zip(X,Y)])
    movement_numerator = log_lik_mov + np.log(prior_a_mov) + np.log(prior_b_mov) + np.log(prior_s_mov)
  
    ratio = np.exp(movement_numerator - current_numerator)
    event = random.uniform(0,1)
    if event <= ratio:
      curr_b = mov_b
      curr_a = mov_a
      current_numerator = movement_numerator
      
  return samples

test2 = lr_mcmc(Y=y,X=x,hops=25_000)
sns.kdeplot([test2[i][0] for i in range(len(test2))],[test2[i][1] for i in range(len(test2))],cmap="inferno",shade=True)  

My plot when running successfully w/o code block change. x-axis = slope, y-axis = y-intercept.

And the error when I change the code

/usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:57: RuntimeWarning: invalid value encountered in double_scalars
/usr/local/lib/python3.6/dist-packages/statsmodels/nonparametric/kernels.py:128: RuntimeWarning: divide by zero encountered in true_divide
  return (1. / np.sqrt(2 * np.pi)) * np.exp(-(Xi - x)**2 / (h**2 * 2.))
/usr/local/lib/python3.6/dist-packages/statsmodels/nonparametric/kernels.py:128: RuntimeWarning: invalid value encountered in true_divide
  return (1. / np.sqrt(2 * np.pi)) * np.exp(-(Xi - x)**2 / (h**2 * 2.))
/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:1483: UserWarning: Warning: converting a masked element to nan.
  self.zmax = float(z.max())
/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:1484: UserWarning: Warning: converting a masked element to nan.
  self.zmin = float(z.min())
/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:1132: RuntimeWarning: invalid value encountered in less
  under = np.nonzero(lev < self.zmin)[0]
/usr/local/lib/python3.6/dist-packages/matplotlib/contour.py:1134: RuntimeWarning: invalid value encountered in greater
  over = np.nonzero(lev > self.zmax)[0]
<matplotlib.axes._subplots.AxesSubplot at 0x7f614fe62ba8>

And when I look at the samples, it's just one b,a,s combination for all 25,00 ierations.

[(1, 1, 1),
 (1, 1, 1),
 (1, 1, 1),
 (1, 1, 1),
 (1, 1, 1),
 (1, 1, 1),
 (1, 1, 1),
 ...
]

Answer 1

There are a few minor issues with your sampler. First, you should always compute the PDF or PMF on a log scale, as for any non-trivial problem you will likely run into overflow/underflow issues (speaking from experience, this happens very quickly). I've switched these functions to a log scale but it's worth checking that I've done this correctly.

The reason your sampler didn't work as you posted it is mainly because you didn't assign curr_s to mov_s (ie, you didn't update the current value if the Metropolis-Hastings ratio was above the uniform value).

For a symmetric non-negative proposal, I've seen people use a log-normal distribution. However you could also adapt the acceptance ratio to be a Metropolis-Hastings ratio rather than a Metropolis ratio (ie, include the ratio of the densities of the proposal as well as the ratio of the posterior, see this for a quick guide). I don't think a truncated uniform distribution is symmetric so the behaviour of your sampler at values of sigma < 0.25 may be incorrect.

You also didn't define your data (X, y) or import all the libraries you used, which is always nice as it makes it easier to load and debug your code.

If I can be so bold as to suggest some simple next steps - it would be pretty straightforward to allow an arbitrary design matrix rather than simply slope and intercept.

As for your question,

What are the benefits of sampling different sigma values?

For the model, $y_i \sim N(X_i\beta, \sigma^2)$, it is natural in most cases to assume that you do not know the magnitude of noise around the conditional mean. Learning this magnitude will enable you to make better predictions using the posterior estimates. This magnitude is controlled by $\sigma^2$. I struggle to think of a situation where the regression coefficients ($\beta$) are unknown, but the magnitude of the residuals $\epsilon_i = y_i - X_i\beta$ is known. Furthermore, to quantify uncertainty around predictions, it's necessary to have a good estimate of the variance around the conditional mean.

import numpy as np
import random
import seaborn as sns

def normalPDF(x,mu,sigma):
  num = np.exp(-1/2*((x-mu)/sigma)**2)
  den = np.sqrt(2*np.pi)*sigma
  return num/den

def invGamma(x,a,b):
  non_zero = int(x>=0)
  func = x**(a-1)*np.exp(-x/b)
  return non_zero*func

def lr_mcmc(X,Y,hops=10_000):
  samples = []
  curr_a = random.gauss(1,1)
  curr_b = random.gauss(2,1)
  curr_s = random.uniform(3,1)

  prior_a_curr = normalPDF(x=curr_a,mu=1,sigma=1)
  prior_b_curr = normalPDF(x=curr_b,mu=2,sigma=1)
  prior_s_curr = invGamma(x=curr_s,a=3,b=1)
  
  log_lik_curr = sum([np.log(normalPDF(x=curr_b*x + curr_a,mu=y,sigma=curr_s)) for x,y in zip(X,Y)])
  current_numerator =  log_lik_curr + np.log(prior_a_curr) + np.log(prior_b_curr) + np.log(prior_s_curr)

  count = 0
  for i in range(hops):
    samples.append((curr_b,curr_a,curr_s))    

    if count == 0: #propose movement to b
      mov_a = curr_a
      mov_b = curr_b + random.uniform(-0.25,0.25) 
      mov_s = curr_s
      count += 1

    elif count == 1: #propose movement to a
      mov_a = curr_a + random.uniform(-0.25,0.25)
      mov_b = curr_b
      mov_s = curr_s
      count += 1

    else: #propose movement to s
      mov_a = curr_a
      mov_b = curr_b
      mov_s = curr_s + random.uniform(-0.25,0.25)
      count = 0

    prior_b_mov = normalPDF(x=mov_b,mu=2,sigma=1)
    prior_a_mov = normalPDF(x=mov_a,mu=1,sigma=1)
    prior_s_mov = invGamma(x=mov_s,a=3,b=1)
    if prior_s_mov <=0: 
      continue #automatically reject because variance cannot equal 0.
    
    log_lik_mov = sum([np.log(normalPDF(x=mov_b*x + mov_a,mu=y,sigma=mov_s)) for x,y in zip(X,Y)])
    movement_numerator = log_lik_mov + np.log(prior_a_mov) + np.log(prior_b_mov) + np.log(prior_s_mov)
  
    ratio = np.exp(movement_numerator - current_numerator)
    event = random.uniform(0,1)
    if event <= ratio:
      curr_b = mov_b
      curr_a = mov_a
      curr_s = mov_s
      current_numerator = movement_numerator
      
  return samples


beta = np.random.normal(0, 1, [1, ])
X = np.random.normal(0, 1, [20, 1])
y = np.matmul(X, beta)

test2 = lr_mcmc(X=X, y=y, hops=25_000)



sns.kdeplot([test2[i][0] for i in range(len(test2))],[test2[i][1] for i in range(len(test2))],cmap="inferno",shade=True)  

plt.show()