MCMC giving strange results for logistic regression example

Hey there, I’m completely new to probabilistic programming, so it may very well be that I made an obvious mistake.

Now to my problem. After reading through the tutorial, Pyro forum and some of my lecture notes, I wanted to implement a logistic regression model based on the bayesian regression tutorial. I finally ended up with the following:

class LogisticRegressionModel(nn.Module):
    def __init__(self, n_features):
        super().__init__()
        self.linear = nn.Linear(n_features, 1)
        self.sigmoid = nn.Sigmoid()

    def forward(self, x):
        return self.sigmoid(self.linear(x))

def model(x, y):
    n_features = x.shape[-1]
    n_samples = x.shape[-2]
    
    loc, scale = torch.zeros(1, n_features),  torch.ones(1, n_features)
    bias_loc, bias_scale = torch.zeros(1), 50 * torch.ones(1)
    w_prior = Normal(loc, scale)
    b_prior = Normal(bias_loc, bias_scale)
    priors = {'linear.weight': w_prior, 'linear.bias': b_prior}
    
    lifted_module = pyro.random_module("module", regression_model, priors)
    
    lifted_reg_model = lifted_module()
    with pyro.plate("map", n_samples):
        model_logits = lifted_reg_model(x).squeeze(-1)
        pyro.sample("obs", Bernoulli(logits=model_logits, validate_args=True), obs=y)

When I now estimate the latent variables through SVI:

optim = Adam({"lr": 0.01})
svi = SVI(model, guide, optim, loss=Trace_ELBO())
regression_model = LogisticRegressionModel(X_train_t.shape[-1])

pyro.clear_param_store()
num_iterations = 10000
for j in range(num_iterations):
    loss = svi.step(X_train_t, y_train_t)
    if j % (num_iterations / 10) == 0:
        print("[iteration %04d] loss: %.4f" % (j + 1, loss / float(X_train_t.shape[-2])))

With a guide that I found on forum here and an automatically generated guide, I get reasonable results, however still considerably worse, when compared to a regular logistic regression model:

def guide(x, y):
    n_features = 29
    
    w_loc = torch.zeros(1, n_features)
    w_log_sig = -3.0 * torch.ones(1, n_features) + 0.05 * torch.randn(1, n_features)
    b_loc = torch.zeros(1)
    b_log_sig = -3.0 * torch.ones(1) + 0.05 * torch.randn(1)
    
    mw_param = pyro.param("guide_mean_weight", w_loc)
    sw_param = softplus(pyro.param("guide_log_scale_weight", w_log_sig))
    mb_param = pyro.param("guide_mean_bias", b_loc)
    sb_param = softplus(pyro.param("guide_log_scale_bias", b_log_sig))
    
    w_dist = Normal(mw_param, sw_param).independent(1)
    b_dist = Normal(mb_param, sb_param).independent(1)
    dists = {'linear.weight': w_dist, 'linear.bias': b_dist}
    
    lifted_module = pyro.random_module("module", regression_model, dists)
    
    return lifted_module()

from pyro.infer.autoguide import AutoDiagonalNormal
guide = AutoDiagonalNormal(model)

y=[guide(None, None)(X_test_t).data.numpy().reshape(-1) for i in range(1000)]
y_pred = np.mean(y,axis=0)
roc_auc_score(y_true=y_test, y_score=y_pred)

I then tried use MCMC sampling using the NUTS kernel:

init_params, potential_fn, transforms, _ = initialize_model(model, model_args=(X_train_t, y_train_t),
                                                                num_chains=1)
nuts_kernel = NUTS(potential_fn=potential_fn)
mcmc = MCMC(nuts_kernel,
            num_samples=5000,
            warmup_steps=500,
            num_chains=1,
            initial_params=init_params,
            transforms=transforms)
mcmc.run(X_train_t, y_train_t)

y=predictive(model, mcmc.get_samples(1), X_test_t, None)
y_pred = y["obs"].mean(axis=0).numpy()

roc_auc_score(y_true=y_test, y_score=y_pred)

When with SVI I got receiver operating characteristic AUC values around .75, the MCMC sampling leads to estimates around .5, thus random. A simple sklearn l2 penalized logistic regression on the same data reaches .84. Now first, there must be a problem somewhere in the MCMC sampling example but I have no clue what it is, could anyone point out the mistake? Second, is it to be expected that a bayesian model has considerably worse performance than its regular frequentist counterpart?

As a final note, a small bug I encountered when initializing multiple chains, it throws the following exception when run in jupyter notebook. I suppose it somehow fails to display the progress bar in this case:

Exception in thread Thread-5:
Traceback (most recent call last):
  File "/home/test/anaconda3/lib/python3.7/threading.py", line 917, in _bootstrap_inner
    self.run()
  File "/home/test/anaconda3/lib/python3.7/threading.py", line 865, in run
    self._target(*self._args, **self._kwargs)
  File "/home/test/anaconda3/lib/python3.7/site-packages/pyro/infer/mcmc/api.py", line 34, in logger_thread
    progress_bars = ProgressBar(warmup_steps, num_samples, disable=disable_progbar, num_bars=num_chains)
  File "/home/test/anaconda3/lib/python3.7/site-packages/pyro/infer/mcmc/logger.py", line 69, in __init__
    position=i, file=sys.stderr, disable=disable)
  File "/home/test/anaconda3/lib/python3.7/site-packages/tqdm/notebook.py", line 201, in __init__
    kwargs['bar_format'] = kwargs['bar_format'].replace('{bar}', '<bar/>')
AttributeError: 'NoneType' object has no attribute 'replace'

Exception ignored in: <function tqdm.__del__ at 0x7fb11d94d598>
Traceback (most recent call last):
  File "/home/test/anaconda3/lib/python3.7/site-packages/tqdm/std.py", line 1039, in __del__
    self.close()
  File "/home/test/anaconda3/lib/python3.7/site-packages/tqdm/notebook.py", line 240, in close
    super(tqdm_notebook, self).close(*args, **kwargs)
  File "/home/test/anaconda3/lib/python3.7/site-packages/tqdm/std.py", line 1215, in close
    if self.disable:
AttributeError: 'tqdm_notebook' object has no attribute 'disable'

Versions:
pyro v: 0.4.1
pytorch v: 1.3.0
Python 3.7
Ubuntu 18.04 LTS

@Mew If you use Bernoulli(logits=model_logits) then it is better to skip sigmoid layer of your LogisticRegressionModel. Otherwise, your LogisticRegressionModel will always give you probs > 0.5. I think that’s why you get a bad estimation of around 0.5. Please let me know if you still get the same problem after removing that layer.

Regarding the progress bar, I guess it is tqdm issue. I am using tqdm=4.32.1 and it works fine in my system. Could you provide your tqdm version so I can take a test for it? In case you are using the latest version of tqdm and still getting progress bar error, could you open an issue at https://github.com/pyro-ppl/pyro/issues? Thanks!

Thanks a lot for the quick answer. You are obviously right with the logits, I should’ve thought about it instead of just copying ;), now results look much better and MCMC is equal to SVI, still not quite like the regular LR but I guess that can be explained by suboptimal priors.

Regarding tqdm I have the most recent version 4.36.1, I’ll quickly check your version to make sure it’s not some other dependency problem and open an issue if it isn’t.