Python emcee.EnsembleSampler() Examples

def sample_representer_points(self):
        self.sampling_acquisition.update(self.model)

        for i in range(5):
            restarts = np.zeros((self.Nb, self.D))
            restarts[0:self.Nb, ] = self.lower + (self.upper - self.lower) \
                                                 * self.rng.uniform(size=(self.Nb, self.D))
            sampler = emcee.EnsembleSampler(
                self.Nb, self.D, self.sampling_acquisition_wrapper)
            # zb are the representer points and lmb are their log EI values
            self.zb, self.lmb, _ = sampler.run_mcmc(restarts, 50)
            if not np.any(np.isinf(self.lmb)):
                break
            else:
                print("Infinity")

        if len(self.zb.shape) == 1:
            self.zb = self.zb[:, None]
        if len(self.lmb.shape) == 1:
            self.lmb = self.lmb[:, None] 

def sample_representer_points(self):
        self.sampling_acquisition.update(self.model)

        start_points = init_random_uniform(self.lower, self.upper, self.Nb)

        sampler = emcee.EnsembleSampler(self.Nb,
                                        self.D,
                                        self.sampling_acquisition_wrapper)

        # zb are the representer points and lmb are their log EI values
        self.zb, self.lmb, _ = sampler.run_mcmc(start_points, 200)

        if len(self.zb.shape) == 1:
            self.zb = self.zb[:, None]
        if len(self.lmb.shape) == 1:
            self.lmb = self.lmb[:, None] 

def __init__(self, model, nwalkers, pool=None,
                 model_call=None):
        try:
            import emcee
        except ImportError:
            raise ImportError("emcee is not installed.")

        if model_call is None:
            model_call = model

        ndim = len(model.variable_params)
        sampler = emcee.EnsembleSampler(nwalkers, ndim,
                                        model_call,
                                        pool=pool)
        # emcee uses it's own internal random number generator; we'll set it
        # to have the same state as the numpy generator
        rstate = numpy.random.get_state()
        sampler.random_state = rstate
        # initialize
        super(EmceeEnsembleSampler, self).__init__(
              sampler, model)
        self._nwalkers = nwalkers 

def fit_hyperparams(self,printOut=False):
    if self.modelp.fitType=='mle':
      spo.minimize(self.neg_log_like, self.model.get_parameter_vector(),
        jac=True)
    elif self.modelp.fitType=='bayes':
      self.nburnin = 200
      nsamp = 200
      nwalkers = 36
      gpdim = len(self.model)
      self.sampler = emcee.EnsembleSampler(nwalkers, gpdim, self.log_post)
      p0 = self.model.get_parameter_vector() + 1e-4*np.random.randn(nwalkers,
        gpdim)
      print 'Running burn-in.'
      p0, _, _ = self.sampler.run_mcmc(p0, self.nburnin)
      print 'Running main chain.'
      self.sampler.run_mcmc(p0, nsamp)
    if printOut:
      print 'Final GP hyperparam (in opt or MCMC chain):'
      print self.model.get_parameter_dict() 

def sample_representer_points(self):
        # Sample representer points only in the
        # configuration space by setting all environmental
        # variables to 1
        D = np.where(self.is_env == 0)[0].shape[0]

        lower = self.lower[np.where(self.is_env == 0)]
        upper = self.upper[np.where(self.is_env == 0)]

        self.sampling_acquisition.update(self.model)

        for i in range(5):
            restarts = np.random.uniform(low=lower,
                                         high=upper,
                                         size=(self.Nb, D))
            sampler = emcee.EnsembleSampler(self.Nb, D,
                                        self.sampling_acquisition_wrapper)

            self.zb, self.lmb, _ = sampler.run_mcmc(restarts, 50)
            if not np.any(np.isinf(self.lmb)):
                break
            else:
                print("Infinity")
        if np.any(np.isinf(self.lmb)):
            raise ValueError("Could not sample valid representer points! LogEI is -infinity")
        if len(self.zb.shape) == 1:
            self.zb = self.zb[:, None]
        if len(self.lmb.shape) == 1:
            self.lmb = self.lmb[:, None]

        # Project representer points to subspace
        proj = np.ones([self.zb.shape[0],
                    self.upper[self.is_env == 1].shape[0]])
        proj *= self.upper[self.is_env == 1].shape[0]
        self.zb = np.concatenate((self.zb, proj), axis=1) 

def emcee_schools_model(data, draws, chains):
    """Schools model in emcee."""
    import emcee

    chains = 10 * chains  # emcee is sad with too few walkers
    y = data["y"]
    sigma = data["sigma"]
    J = data["J"]  # pylint: disable=invalid-name
    ndim = J + 2

    pos = np.random.normal(size=(chains, ndim))
    pos[:, 1] = np.absolute(pos[:, 1])  #  pylint: disable=unsupported-assignment-operation

    if emcee_version() < 3:
        sampler = emcee.EnsembleSampler(chains, ndim, _emcee_lnprob, args=(y, sigma))
        # pylint: enable=unexpected-keyword-arg
        sampler.run_mcmc(pos, draws)
    else:
        here = os.path.dirname(os.path.abspath(__file__))
        data_directory = os.path.join(here, "saved_models")
        filepath = os.path.join(data_directory, "reader_testfile.h5")
        backend = emcee.backends.HDFBackend(filepath)  # pylint: disable=no-member
        backend.reset(chains, ndim)
        # pylint: disable=unexpected-keyword-arg
        sampler = emcee.EnsembleSampler(
            chains, ndim, _emcee_lnprob, args=(y, sigma), backend=backend
        )
        # pylint: enable=unexpected-keyword-arg
        sampler.run_mcmc(pos, draws, store=True)
    return sampler


# pylint:disable=no-member,no-value-for-parameter,invalid-name 

def get_samples(self, n_samples, log_p_function, burn_in_steps=50, n_steps=100):
        """
        Generates samples.

        Parameters:
            n_samples - number of samples to generate
            log_p_function - a function that returns log density for a specific sample
            burn_in_steps - number of burn-in steps for sampling

        Returns a tuple of two array: (samples, log_p_function values for samples)
        """
        X_init = self.space.sample_uniform(n_samples)
        sampler = emcee.EnsembleSampler(n_samples, X_init.shape[1], log_p_function)

        # Burn-In
        state = list(sampler.run_mcmc(X_init, burn_in_steps)) # compatible with both emcee 2 and 3
        samples = state[0]
        samples_log = state[1]

        # MCMC Sampling
        state = list(sampler.run_mcmc(samples, n_steps))
        samples = state[0]
        samples_log = state[1]

        # make sure we have an array of shape (n samples, space input dim)
        if len(samples.shape) == 1:
            samples = samples.reshape(-1, 1)
        samples_log = samples_log.reshape(-1, 1)

        return samples, samples_log 

def __init__(self, ntemps, nwalkers, dim, logl, logp, threads=1,
                 pool=None, betas=None):
        self.logl = logl
        self.logp = logp

        self.ntemps = ntemps
        self.nwalkers = nwalkers
        self.dim = dim

        self._chain = None
        self._lnprob = None
        self._lnlikelihood = None

        if betas is None:
            self._betas = self.exponential_beta_ladder(ntemps)
        else:
            self._betas = betas

        self.nswap = np.zeros(ntemps, dtype=np.float)
        self.nswap_accepted = np.zeros(ntemps, dtype=np.float)

        self.pool = pool
        if threads > 1 and pool is None:
            self.pool = multi.Pool(threads)

        self.samplers = [em.EnsembleSampler(nwalkers, dim,
                                            PTPost(logl, logp, b),
                                            pool=self.pool)
                                    for b in self.betas] 

def _initialise_sampler(self):
        self._sampler = self.emcee.EnsembleSampler(**self.sampler_init_kwargs)
        self._init_chain_file() 

def __init__(self, model, nwalkers,
                 checkpoint_interval=None, checkpoint_signal=None,
                 logpost_function=None, nprocesses=1, use_mpi=False):

        self.model = model
        # create a wrapper for calling the model
        if logpost_function is None:
            logpost_function = 'logposterior'
        model_call = models.CallModel(model, logpost_function)

        # Set up the pool
        if nprocesses > 1:
            # these are used to help paralleize over multiple cores / MPI
            models._global_instance = model_call
            model_call = models._call_global_model
        pool = choose_pool(mpi=use_mpi, processes=nprocesses)
        if pool is not None:
            pool.count = nprocesses

        # set up emcee
        self.nwalkers = nwalkers
        ndim = len(model.variable_params)
        self._sampler = emcee.EnsembleSampler(nwalkers, ndim, model_call,
                                              pool=pool)
        # emcee uses it's own internal random number generator; we'll set it
        # to have the same state as the numpy generator
        rstate = numpy.random.get_state()
        self._sampler.random_state = rstate
        self._checkpoint_interval = checkpoint_interval
        self._checkpoint_signal = checkpoint_signal