0
$\begingroup$

Trying to understand the pyro example here: https://pyro.ai/examples/svi_part_i.html which starts with a Beta(10,10) prior, adds 10 Bernoulli likelihood datapoints with a 6,4 split. The analytic solution gives a posterior of Beta(16,14).

The guide starts at Beta(15,15) and after 2000 steps and 11 seconds of processing I get an estimated posterior of Beta(16.10,13.76) where I expected it to lock onto the exact answer. I also tried starting with the correct answer for the guide and processing for 20,000 steps taking 110 seconds.

\begin{array} \ Guide & steps & Posterior \\ Beta(15,15) & 2000 & Beta(16.10,13.76) \\ Beta(15,15) & 2000 & Beta(15.94,14.10) \\ Beta(15,15) & 2000 & Beta(16.05,13.83) \\ Beta(15,15) & 20000 & Beta(16.13,13.80) \\ Beta(16,14) & 2000 & Beta(16.04,13.97) \\ Beta(16,14) & 2000 & Beta(15.99,14.02) \\ Beta(16,14) & 20000 & Beta(15.82,14.09) \\ \end{array} This is a toy problem and yet there are several issues:

  • Unable to identify optimal solution
  • More steps don't improve results

Are my expectations too high? Why is it unable to stick with the optimal solution?

from __future__ import print_function
import math
import os
import torch
import torch.distributions.constraints as constraints
import pyro
from pyro.optim import Adam
from pyro.infer import SVI, Trace_ELBO
import pyro.distributions as dist

# this is for running the notebook in our testing framework
smoke_test = ('CI' in os.environ)
n_steps = 2 if smoke_test else 2000

# enable validation (e.g. validate parameters of distributions)
assert pyro.__version__.startswith('0.3.1')
pyro.enable_validation(True)

# clear the param store in case we're in a REPL
pyro.clear_param_store()

# create some data with 6 observed heads and 4 observed tails
data = []
for _ in range(6):
    data.append(torch.tensor(1.0))
for _ in range(4):
    data.append(torch.tensor(0.0))

def model(data):
    # define the hyperparameters that control the beta prior
    alpha0 = torch.tensor(10.0)
    beta0 = torch.tensor(10.0)
    # sample f from the beta prior
    f = pyro.sample("latent_fairness", dist.Beta(alpha0, beta0))
    # loop over the observed data
    for i in range(len(data)):
        # observe datapoint i using the bernoulli likelihood
        pyro.sample("obs_{}".format(i), dist.Bernoulli(f), obs=data[i])

def guide(data):
    # register the two variational parameters with Pyro
    # - both parameters will have initial value 15.0.
    # - because we invoke constraints.positive, the optimizer
    # will take gradients on the unconstrained parameters
    # (which are related to the constrained parameters by a log)
    alpha_q = pyro.param("alpha_q", torch.tensor(15.0),
                         constraint=constraints.positive)
    beta_q = pyro.param("beta_q", torch.tensor(15.0),
                        constraint=constraints.positive)
    # sample latent_fairness from the distribution Beta(alpha_q, beta_q)
    pyro.sample("latent_fairness", dist.Beta(alpha_q, beta_q))

# setup the optimizer
adam_params = {"lr": 0.0005, "betas": (0.90, 0.999)}
optimizer = Adam(adam_params)

# setup the inference algorithm
svi = SVI(model, guide, optimizer, loss=Trace_ELBO())

# do gradient steps
for step in range(n_steps):
    svi.step(data)
    if step % 100 == 0:
        print('.', end='')

# grab the learned variational parameters
alpha_q = pyro.param("alpha_q").item()
beta_q = pyro.param("beta_q").item()

# here we use some facts about the beta distribution
# compute the inferred mean of the coin's fairness
inferred_mean = alpha_q / (alpha_q + beta_q)
# compute inferred standard deviation
factor = beta_q / (alpha_q * (1.0 + alpha_q + beta_q))
inferred_std = inferred_mean * math.sqrt(factor)

print("\nbased on the data and our prior belief, the fairness " +
      "of the coin is %.3f +- %.3f" % (inferred_mean, inferred_std))

I lack the reputation to create a pyro tag.

$\endgroup$
  • 1
    $\begingroup$ Have you tried more iterations? $\endgroup$ – Tim Apr 19 at 15:13
  • $\begingroup$ 20000 steps moved the result further from the truth than 2000 as shown in the table. $\endgroup$ – Andrew Apr 20 at 4:07
3
$\begingroup$

As always when using SGD, care has to be taken with the learning rate. By decaying learning rate from 1E-3 to 8E-5 over 5000 iterations ($\gamma = 0.995$) using the SGD optimizer with momentum 0.9, I was able to get these results from 3 runs of the algorithm: for $\alpha$ and $\beta$ respectively:

15.997169494628906 13.988293647766113
15.993146896362305 14.00241756439209
15.991819381713867 14.001907348632812

I also set the batch size to 4096 to reduce the the variance of the gradient. Since most of the runtime is in the library overhead, this barely increases the runtime over the default batch size of 1.

Iteration count divided by 10 on the x-axis, $\beta_q$ on the $y$ axis. As the learning rate is still decaying, more iterations would continue to improve the convergence.

$\endgroup$

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.