I'm trying to reproduce Figure 2 from this paper. In summary, I have the regression model

$$Y \sim N(\mu, \sigma^2) \\ \mu = \alpha + \beta_1X_1 + \beta_2X_2 $$

The prior distributions for the parameters are

$$ \alpha, \beta_1 , \beta_2 \sim N(0, 10)\\ \sigma \sim \vert N(0, 1)\vert $$

The goal is to compute posterior distributions for parameters based on the (synthetically generated) data using vanilla MCMC (in contrast to the NUTS sampler as done in the paper).

I'm defining 'state' $s$ as the tuple $s=(\alpha, \beta_1, \beta_2, \sigma)$. The log-likelihood of the data given the parameters is

$$ \log P(Y \vert s) =-\frac{M}{2}\log(2\pi\sigma^2) - \frac{1}{2\sigma^2}\sum_{i=1}^M\left(y_i -w^Tx_i\right)^2 $$

where $w = [\alpha, \beta_1, \beta_2]$ is the parameter vector, $x_i$ is the $i$th row of the (synthetic) data, with the bias term. I.e. $x_i = [1, x_{1i}, x_{2i}]$. In the code below, you'll see the first term computed as $-M\log\sigma$ since the constants will cancel out in the eventual subtraction in the acceptance calculation.

The log-prior is simply the sum of the logarithms of the PDFs of the parameters as described above (computed in the code below).

I'm using the random-walk MCMC where the proposed next state is chosen with the current state as the mean and a fixed variance vector as the scale hyper parameter.

$$ s^{t+1} \sim N(s^t, \nu) $$

with $\nu = [0.07, 0.25, 0.05, 0.02]$. With these manually tweaked, current values, I'm getting incorrect results. The results are problematic in two respects: (1) They converge to wrong values and (2) I'm generating correlated samples. I've experimented with the $\nu$ parameter for the whole weekend and have not had great success getting it to work.


As can be seen all of my trace plots show severe sample correlation problem. The posterior modes also converge to incorrect values (true values should be $[\alpha=1, \beta_1=1, \beta_2=2.5, \sigma=1]$.

I feel like I'm missing something obvious since this 4-dimensional problem should be solvable using the Metropolis algorithm. I'm including the full code below. Any help is appreciated.

import numpy as np
from scipy.stats import norm, halfnorm
import matplotlib.pyplot as pl
from tqdm import tqdm

def trace():
    import ipdb; ipdb.set_trace()

def loglik(state, X, Y):
    M = X.shape[0]
    coefs, sig = state[:-1], state[-1]
    err = Y - np.dot(X, coefs)   
    term1 = -M * np.log(sig)
    term2 = -np.sum(err * err) / (2.0 * sig * sig)    
    return term1 + term2    

def proposalfunc(state):    
    return np.random.normal(loc=state, scale=[0.07, 0.25, 0.05, 0.02])        

def logprior(state, distributions):
    return sum(np.log(D.pdf(x)) for x, D in zip(state, distributions))

def generate_data(size):
    alpha, beta1, beta2, sig = 1.0, 1.0, 2.5, 1.0    
    X1 = np.linspace(0, 1.0, size)
    X2 = np.linspace(0, 0.2, size)
    ONES = np.ones(X1.shape)   
    noise = np.random.normal(loc=0, scale=sig, size=size)
    X = np.vstack((ONES, X1, X2)).T
    Y = np.dot(X, [alpha, beta1, beta2]) + noise   
    return X, Y

def plot(chain, skip):
    param_names = ['alpha', 'beta1', 'beta2', 'sigma']
    fig, axes = pl.subplots(4, 2, figsize=(10, 6))
    histargs = dict(normed=True, histtype='step', bins=50, color='r', alpha=0.7)
    for i, pname in zip(range(chain.shape[1]), param_names):
        param = chain[skip:, i]
        axes[i, 0].hist(param, **histargs)        
        axes[i, 1].plot(param, alpha=0.5)        
        axes[i, 0].set_title(pname)

def mcmc_regression():    
    np.random.seed(123) # Uncomment to reproduce the plot exactly
    X, Y = generate_data(100)    
    alpha_dist = beta1_dist = beta2_dist = norm(0, 10)
    sigma_dist = halfnorm(0, 1)
    dists = (alpha_dist, beta1_dist, beta2_dist, sigma_dist)

    nburnin = 10000
    niter = 5000 + nburnin
    ncomps = 4
    chain = np.zeros(shape=(niter, ncomps))
    chain[0, :] = np.random.normal(size=ncomps)

    for i in tqdm(range(niter - 1)):                
        v = chain[i]
        log_posterior_old = loglik(v, X, Y) + logprior(v, dists)         
        proposal = proposalfunc(v)
        log_posterior_new = loglik(proposal, X, Y) + logprior(proposal, dists)        
        a = min(0.0, log_posterior_new - log_posterior_old)
        if np.random.random() < np.exp(a):            
            chain[i + 1, :] = proposal
            chain[i + 1, :] = v

    plot(chain, nburnin)

if __name__ == '__main__':

EDIT: I forgot to mention a few things in the OP.

  1. I'm able to get a univariate Bayesian MCMC linear regression to work. But am a bit lost on this bivariate case.

  2. The results are quite sensitive to the randomizer seed (which further indicates that something is wrong at a very basic level).

  3. I'm aware that HMC and NUTS can help mitigate the sample correlation problem, but I want to get a good understanding of the simple case first.

  • 1
    $\begingroup$ I may have found one error. For $\sigma$ only positive values are allowed, but the proposal draws from a normal distribution, which may propose negative values. Such proposal values should immediately be rejected. $\endgroup$ – Greenparker Apr 23 '18 at 6:23
  • 1
    $\begingroup$ Also, what is the acceptance rate for your chain? Ideally $\nu$ should be chosen so that the acceptance rate is close to .234. $\endgroup$ – Greenparker Apr 23 '18 at 6:24
  • 1
    $\begingroup$ Also your initialization is from a normal distribution, which may also take negative values for $\sigma$. $\endgroup$ – Greenparker Apr 23 '18 at 6:40
  • $\begingroup$ I am curious; was any of this helpful? $\endgroup$ – Greenparker Apr 24 '18 at 23:28
  • $\begingroup$ @Greenparker Thank you for your comments and sincere apologies for not updating sooner. I've only implemented halfnormal for $\sigma$ initialization. It is helped only marginally. I think the acceptance rate could be crucial. I'll implement that next. $\endgroup$ – RDK May 22 '18 at 14:20

Your Answer

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

Browse other questions tagged or ask your own question.