# Simple MCMC Linear Regression

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

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.

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`

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