code: import numpy as np
  
  from ..arraystep import Competence
  from pymc3_ext.vartypes import discrete_types
  from pymc3_ext.step_methods.hmc.integration import IntegrationError
  from pymc3_ext.step_methods.hmc.base_hmc import BaseHMC, HMCStepData, DivergenceInfo
  
  
  __all__ = ['HamiltonianMC']
  
  
  def unif(step_size, elow=.85, ehigh=1.15):
      return np.random.uniform(elow, ehigh) * step_size
  
  
  class HamiltonianMC(BaseHMC):
      R"""A sampler for continuous variables based on Hamiltonian mechanics.
  
      See NUTS sampler for automatically tuned stopping time and step size scaling.
      """
  
      name = 'hmc'
      default_blocked = True
      generates_stats = True
      stats_dtypes = [{
          'step_size': np.float64,
          'n_steps': np.int64,
          'tune': np.bool,
          'step_size_bar': np.float64,
          'accept': np.float64,
          'diverging': np.bool,
          'energy_error': np.float64,
          'energy': np.float64,
          'max_energy_error': np.float64,
          'path_length': np.float64,
          'accepted': np.bool,
          'model_logp': np.float64,
      }]
  
      def __init__(self, vars=None, path_length=2., **kwargs):
          """Set up the Hamiltonian Monte Carlo sampler.
  
          Parameters
          ----------
          vars : list of theano variables
          path_length : float, default=2
              total length to travel
          step_rand : function float -> float, default=unif
              A function which takes the step size and returns an new one used to
              randomize the step size at each iteration.
          step_scale : float, default=0.25
              Initial size of steps to take, automatically scaled down
              by 1/n**(1/4).
          scaling : array_like, ndim = {1,2}
              The inverse mass, or precision matrix. One dimensional arrays are
              interpreted as diagonal matrices. If `is_cov` is set to True,
              this will be interpreded as the mass or covariance matrix.
          is_cov : bool, default=False
              Treat the scaling as mass or covariance matrix.
          potential : Potential, optional
              An object that represents the Hamiltonian with methods `velocity`,
              `energy`, and `random` methods. It can be specified instead
              of the scaling matrix.
          target_accept : float, default 0.65
              Adapt the step size such that the average acceptance
              probability across the trajectories are close to target_accept.
              Higher values for target_accept lead to smaller step sizes.
              Setting this to higher values like 0.9 or 0.99 can help
              with sampling from difficult posteriors. Valid values are
              between 0 and 1 (exclusive). Default of 0.65 is from (Beskos et.
              al. 2010, Neal 2011). See Hoffman and Gelman's "The No-U-Turn
              Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte
              Carlo" section 3.2 for details.
          gamma : float, default .05
          k : float, default .75
              Parameter for dual averaging for step size adaptation. Values
              between 0.5 and 1 (exclusive) are admissible. Higher values
              correspond to slower adaptation.
          t0 : float > 0, default 10
              Parameter for dual averaging. Higher values slow initial
              adaptation.
          adapt_step_size : bool, default=True
              Whether step size adaptation should be enabled. If this is
              disabled, `k`, `t0`, `gamma` and `target_accept` are ignored.
          model : pymc3.Model
              The model
          **kwargs : passed to BaseHMC
          """
          kwargs.setdefault('step_rand', unif)
          kwargs.setdefault('target_accept', 0.65)
          super().__init__(vars, **kwargs)
          self.path_length = path_length
  
      def _hamiltonian_step(self, start, p0, step_size):
          n_steps = max(1, int(self.path_length / step_size))
  
          energy_change = -np.inf
          state = start
          div_info = None
          try:
              for _ in range(n_steps):
                  state = self.integrator.step(step_size, state)
          except IntegrationError as e:
              div_info = DivergenceInfo('Divergence encountered.', e, state)
          else:
              if not np.isfinite(state.energy):
                  div_info = DivergenceInfo(
                      'Divergence encountered, bad energy.', None, state)
              energy_change = start.energy - state.energy
              if np.abs(energy_change) > self.Emax:
                  div_info = DivergenceInfo(
                      'Divergence encountered, large integration error.',
                      None, state)
  
          accept_stat = min(1, np.exp(energy_change))
  
          if div_info is not None or np.random.rand() >= accept_stat:
              end = start
              accepted = False
          else:
              end = state
              accepted = True
  
          stats = {
              'path_length': self.path_length,
              'n_steps': n_steps,
              'accept': accept_stat,
              'energy_error': energy_change,
              'energy': state.energy,
              'accepted': accepted,
              'model_logp': state.model_logp,
          }
          return HMCStepData(end, accept_stat, div_info, stats)
  
      @staticmethod
      def competence(var, has_grad):
          """Check how appropriate this class is for sampling a random variable."""
          if var.dtype in discrete_types or not has_grad:
              return Competence.INCOMPATIBLE
          return Competence.COMPATIBLE
   import pytest
  from numpy import ma, array
  import numpy
  import pandas as pd
  from pymc3_ext import Model, Normal, sample_prior_predictive, sample, ImputationWarning
  
  def test_missing():
      data = ma.masked_values([1, 2, -1, 4, -1], value=-1)
      with Model() as model:
          x = Normal('x', 1, 1)
          with pytest.warns(ImputationWarning):
              Normal('y', x, 1, observed=data)
  
      y_missing, = model.missing_values
      assert y_missing.tag.test_value.shape == (2,)
  
      model.logp(model.test_point)
  
      with model:
          prior_trace = sample_prior_predictive()
      assert set(['x', 'y']) <= set(prior_trace.keys())
  
  
  def test_missing_pandas():
      data = pd.DataFrame([1, 2, numpy.nan, 4, numpy.nan])
      with Model() as model:
          x = Normal('x', 1, 1)
          with pytest.warns(ImputationWarning):
              Normal('y', x, 1, observed=data)
  
      y_missing, = model.missing_values
      assert y_missing.tag.test_value.shape == (2,)
  
      model.logp(model.test_point)
  
      with model:
          prior_trace = sample_prior_predictive()
      assert set(['x', 'y']) <= set(prior_trace.keys())
  
  def test_missing_with_predictors():
      predictors = array([0.5, 1, 0.5, 2, 0.3])
      data = ma.masked_values([1, 2, -1, 4, -1], value=-1)
      with Model() as model:
          x = Normal('x', 1, 1)
          with pytest.warns(ImputationWarning):
              Normal('y', x * predictors, 1, observed=data)
  
      y_missing, = model.missing_values
      assert y_missing.tag.test_value.shape == (2,)
  
      model.logp(model.test_point)
  
      with model:
          prior_trace = sample_prior_predictive()
      assert set(['x', 'y']) <= set(prior_trace.keys())
  
  
  def test_missing_dual_observations():
      with Model() as model:
          obs1 = ma.masked_values([1, 2, -1, 4, -1], value=-1)
          obs2 = ma.masked_values([-1, -1, 6, -1, 8], value=-1)
          beta1 = Normal('beta1', 1, 1)
          beta2 = Normal('beta2', 2, 1)
          latent = Normal('theta', shape=5)
          with pytest.warns(ImputationWarning):
              ovar1 = Normal('o1', mu=beta1 * latent, observed=obs1)
          with pytest.warns(ImputationWarning):
              ovar2 = Normal('o2', mu=beta2 * latent, observed=obs2)
  
          prior_trace = sample_prior_predictive()
          assert set(['beta1', 'beta2', 'theta', 'o1', 'o2']) <= set(prior_trace.keys())
          sample()
  
  def test_internal_missing_observations():
      with Model() as model:
          obs1 = ma.masked_values([1, 2, -1, 4, -1], value=-1)
          obs2 = ma.masked_values([-1, -1, 6, -1, 8], value=-1)
          with pytest.warns(ImputationWarning):
              theta1 = Normal('theta1', mu=2, observed=obs1)
          with pytest.warns(ImputationWarning):
              theta2 = Normal('theta2', mu=theta1, observed=obs2)
  
          prior_trace = sample_prior_predictive()
          assert set(['theta1', 'theta2']) <= set(prior_trace.keys())
          sample()
  
apis: numpy.random.uniform numpy.abs numpy.exp numpy.random.rand numpy.isfinite
apis: numpy.array pandas.DataFrame numpy.ma.masked_values