docs: def variational_lower_bound(params, t, logprob, sampler, log_density,
                              num_samples, rs):
      """Provides a stochastic estimate of the variational lower bound,
         for any variational family and model density."""
      samples = sampler(params, num_samples, rs)
      log_qs = log_density(params, samples)
      log_ps = logprob(samples, t)
      log_ps = np.reshape(log_ps, (num_samples, -1))
      log_qs = np.reshape(log_qs, (num_samples, -1))
      return np.mean(log_ps - log_qs)