7
$\begingroup$

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. graph

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),
 ...
]
$\endgroup$
6
  • 1
    $\begingroup$ Can you explain how the sampler is hindered? Does 1 iteration takes a long time? High autocorelation in chains? Anything else? As a side note you should estimate $\sigma$ anyway since it is a very important part of the model: I cannot produce reliable prediction intervals by conditioning on an arbitrary value of $\sigma$. $\endgroup$
    – jcken
    Sep 14, 2020 at 18:51
  • $\begingroup$ @jcken, Sorry for the confusion. When I propose changes to sigma, concurrently with changes to b and a, the movement is rejected. My last attempt resulted in 25,000 samples of the same point. Perhaps I should adjust variables one at a time...? $\endgroup$
    – jbuddy_13
    Sep 14, 2020 at 21:07
  • $\begingroup$ Yes, changing one at a time can help a lot. Also I noticed you are computing the normal pdf incorrectly! $\endgroup$
    – jcken
    Sep 15, 2020 at 6:02
  • $\begingroup$ @jcken, I omitted the scaling factor in the denominator, but retained the exponentiated term as is. This is very common in MH, as the PDF isn't proposing samples; rather it's returning a likelihood. And for MH purposes, the likelihoods need not be scaled, they just need to be evaluated in a consistent manner! (Long story short, I don't think I made any errors on the exp term.) $\endgroup$
    – jbuddy_13
    Sep 15, 2020 at 16:19
  • 1
    $\begingroup$ I assume the error is just because you're trying to do KDE on a single value. Your sampler has a few issues that I've highlighted/fixed below. $\endgroup$ Sep 15, 2020 at 17:58

1 Answer 1

6
$\begingroup$

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. 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.

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()
```
$\endgroup$
3
  • 1
    $\begingroup$ Awesome, thank you! This honestly helped so much :) $\endgroup$
    – jbuddy_13
    Sep 16, 2020 at 12:09
  • $\begingroup$ I am a beginner in learning mcmc and I found your post really helpful. Could you please post the data you've used in this code? $\endgroup$
    – nancy
    May 12 at 15:56
  • $\begingroup$ @nancy, the data in the code is just simulated data that's generated using the np.random code at the end $\endgroup$ May 12 at 20:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.