Markov chain Monte Carlo (MCMC) methods are considered the gold standard of Bayesian inference; under suitable conditions and in the limit of infinitely many draws they generate samples from the true posterior distribution. HMC (Neal, 2011) uses gradients of the model's log-density function to propose samples, allowing it to exploit posterior geometry. However, it is computationally more expensive than variational inference and relatively sensitive to tuning.

  num_results = 100,
  num_warmup_steps = 50,
  num_leapfrog_steps = 15,
  initial_state = NULL,
  initial_step_size = NULL,
  chain_batch_shape = list(),
  num_variational_steps = 150,
  variational_optimizer = NULL,
  variational_sample_size = 5,
  seed = NULL,
  name = NULL



float tensor of shape concat([sample_shape, model.batch_shape, [num_timesteps, 1]]) where sample_shape corresponds to i.i.d. observations, and the trailing [1] dimension may (optionally) be omitted if num_timesteps > 1. May optionally be an instance of sts_masked_time_series, which includes a mask tensor to specify timesteps with missing observations.


An instance of StructuralTimeSeries representing a time-series model. This represents a joint distribution over time-series and their parameters with batch shape [b1, ..., bN].


Integer number of Markov chain draws. Default value: 100.


Integer number of steps to take before starting to collect results. The warmup steps are also used to adapt the step size towards a target acceptance rate of 0.75. Default value: 50.


Integer number of steps to run the leapfrog integrator for. Total progress per HMC step is roughly proportional to step_size * num_leapfrog_steps. Default value: 15.


Optional Python list of Tensors, one for each model parameter, representing the initial state(s) of the Markov chain(s). These should have shape tf$concat(list(chain_batch_shape, param$prior$batch_shape, param$prior$event_shape)). If NULL, the initial state is set automatically using a sample from a variational posterior. Default value: NULL.


list of tensors, one for each model parameter, representing the step size for the leapfrog integrator. Must broadcast with the shape of initial_state. Larger step sizes lead to faster progress, but too-large step sizes make rejection exponentially more likely. If NULL, the step size is set automatically using the standard deviation of a variational posterior. Default value: NULL.


Batch shape (list or int) of chains to run in parallel. Default value: list() (i.e., a single chain).


int number of steps to run the variational optimization to determine the initial state and step sizes. Default value: 150.


Optional tf$train$Optimizer instance to use in the variational optimization. If NULL, defaults to tf$train$AdamOptimizer(0.1). Default value: NULL.


integer number of Monte Carlo samples to use in estimating the variational divergence. Larger values may stabilize the optimization, but at higher cost per step in time and memory. Default value: 1.


integer to seed the random number generator.


name prefixed to ops created by this function. Default value: NULL (i.e., 'fit_with_hmc').


list of:

  • samples: list of Tensors representing posterior samples of model parameters, with shapes [concat([[num_results], chain_batch_shape, param.prior.batch_shape, param.prior.event_shape]) for param in model.parameters].

  • kernel_results: A (possibly nested) list of Tensors representing internal calculations made within the HMC sampler.


This method attempts to provide a sensible default approach for fitting StructuralTimeSeries models using HMC. It first runs variational inference as a fast posterior approximation, and initializes the HMC sampler from the variational posterior, using the posterior standard deviations to set per-variable step sizes (equivalently, a diagonal mass matrix). During the warmup phase, it adapts the step size to target an acceptance rate of 0.75, which is thought to be in the desirable range for optimal mixing (Betancourt et al., 2014).


See also


# \donttest{ observed_time_series <- rep(c(3.5, 4.1, 4.5, 3.9, 2.4, 2.1, 1.2), 5) + rep(c(1.1, 1.5, 2.4, 3.1, 4.0), each = 7) %>% tensorflow::tf$convert_to_tensor(dtype = tensorflow::tf$float64)
#> Error in py_call_impl(callable, dots$args, dots$keywords): InternalError: CUDA runtime implicit initialization on GPU:0 failed. Status: out of memory #> #> Detailed traceback: #> File "/home/key/.virtualenvs/tf2.2/lib64/python3.7/site-packages/tensorflow/python/framework/", line 1283, in convert_to_tensor_v2 #> as_ref=False) #> File "/home/key/.virtualenvs/tf2.2/lib64/python3.7/site-packages/tensorflow/python/framework/", line 1341, in convert_to_tensor #> ret = conversion_func(value, dtype=dtype, name=name, as_ref=as_ref) #> File "/home/key/.virtualenvs/tf2.2/lib64/python3.7/site-packages/tensorflow/python/framework/", line 321, in _constant_tensor_conversion_function #> return constant(v, dtype=dtype, name=name) #> File "/home/key/.virtualenvs/tf2.2/lib64/python3.7/site-packages/tensorflow/python/framework/", line 262, in constant #> allow_broadcast=True) #> File "/home/key/.virtualenvs/tf2.2/lib64/python3.7/site-packages/tensorflow/python/framework/", line 270, in _constant_impl #> t = convert_to_eager_tensor(value, ctx, dtype) #> File "/home/key/.virtualenvs/tf2.2/lib64/python3.7/site-packages/tensorflow/python/framework/", line 95, in convert_to_eager_tensor #> ctx.ensure_initialized() #> File "/home/key/.virtualenvs/tf2.2/lib64/python3.7/site-packages/tensorflow/python/eager/", line 515, in ensure_initialized #> context_handle = pywrap_tfe.TFE_NewContext(opts)
day_of_week <- observed_time_series %>% sts_seasonal(num_seasons = 7)
#> Error in eval(lhs, parent, parent): object 'observed_time_series' not found
local_linear_trend <- observed_time_series %>% sts_local_linear_trend()
#> Error in eval(lhs, parent, parent): object 'observed_time_series' not found
model <- observed_time_series %>% sts_sum(components = list(day_of_week, local_linear_trend))
#> Error in eval(lhs, parent, parent): object 'observed_time_series' not found
states_and_results <- observed_time_series %>% sts_fit_with_hmc( model, num_results = 10, num_warmup_steps = 5, num_variational_steps = 15)
#> Error in eval(lhs, parent, parent): object 'observed_time_series' not found
# }