Simple MCMC Linear Regression

1

Hello,

When using the example below I get a model specification incorrect - trace log pdf is NaN or Inf, wondering if someone could please tell me what I’m doing wrong here? It would be much appreciated!

import torch
import pyro
from pyro.distributions import Normal
from pyro.infer.mcmc import NUTS, MCMC
from pyro.infer import EmpiricalMarginal

X = torch.randn((100,))
b = torch.tensor(0.25)
noise = torch.randn((100,))
y = X*b + noise

def model(data):
  b_est  = pyro.sample("beta", Normal(loc=torch.zeros(1), scale=torch.ones(1)))
  s2_est = pyro.sample("variance", Normal(loc=torch.zeros(1), scale=torch.ones(1)))
  eps = pyro.sample("noise", Normal(loc=torch.zeros(100), scale=torch.ones(100)))
  y_hat = pyro.sample("yhat", Normal(loc=data*b_est + eps, scale=s2_est))
  return y_hat

nuts_kernel = NUTS(model, step_size=0.09)
mcmc_run = MCMC(nuts_kernel, num_samples=500, warmup_steps=100).run(X)
posterior = EmpiricalMarginal(mcmc_run, "beta")
posterior.mean
2

A few observations:

  • I think you missed including your observed data y in your model, via the obs keyword when you sample y_hat.
  • You should also constrain your “variance” to be positive. The simplest way would be to use a HalfCauchy distribution instead of a Normal.
  • I don’t think you would want to sample eps, since the s2_est factor would account for the noise in your observed data.
3

thanks @neerajprad, that works

How would one get feedback on what parameters to change/tune in the model so to get better performance? For the problem above - though it’s not widely terrible the estimated coefficients range from 0.14 - 0.4, this is with your recommendation of removing eps

4

How are you measuring the range? The mean should be close to 0.25, but you will need more data points to shrink your posterior distribution. This is expected as the half-cauchy prior on the scale parameter is a very weak one.

e.g. code below for 10k data points. You can also get a narrower posterior distribution by making stronger assumptions in your model (which may not always be applicable), for instance, by fixing s2_est = 0.05.

import torch
import pyro
from pyro.distributions import Normal, HalfCauchy
from pyro.infer.mcmc import NUTS, MCMC
from pyro.infer import EmpiricalMarginal
import seaborn as sns
import matplotlib.pyplot as plt


def model(X, y):
    b_est = pyro.sample("beta", Normal(loc=torch.zeros(1), scale=torch.ones(1)))
    s2_est = pyro.sample("variance", HalfCauchy(loc=torch.zeros(1), scale=torch.ones(1)))
    y_hat = pyro.sample("yhat", Normal(loc=b_est * X , scale=s2_est), obs=y)
    return y_hat


X = torch.randn((10000,))
b = torch.tensor(0.25)
noise = torch.randn((10000,))
y = X * b + noise
nuts_kernel = NUTS(model, adapt_step_size=True)
mcmc_run = MCMC(nuts_kernel, num_samples=1000, warmup_steps=400).run(X, y)
posterior = EmpiricalMarginal(mcmc_run, "beta")
print(posterior.mean)
sns.distplot(posterior.get_samples_and_weights()[0])
plt.show()