docs: def infer(self, data, initial_proposal=None, full_output=False,**kwargs):
  
          """
          Infer the model parameters, given the data.
          auto_convergence=True,
          walkers=100, burn=2000, sample=2000, minimum_sample=2000,
          convergence_check_frequency=1000, a=2.0, threads=1,
  
          """
  
          # Apply data masks now so we don't have to do it on the fly.
          data, pixels_affected = self._apply_data_mask(data)
  
          # Any channels / parameters to ignore?
          matched_channels, missing_channels, ignore_parameters \
              = self._match_channels_to_data(data)
          parameters = [p for p in self.parameters if p not in ignore_parameters]
          #parameters = list(set(self.parameters).difference(ignore_parameters))
  
          logger.debug("Inferring {0} parameters: {1}".format(len(parameters),
              ", ".join(parameters)))
  
          # What sampling behaviour will we have?
          # - Auto-convergence:
          #       + Sample for `minimum_sample` (default 2000, 200 walkers)
          #       + Calculate the maximum exponential autocorrelation time for
          #         all parameters
          #       + For the rest of the chain, calculate the autocorrelation time
          #       + Ensure that the number of samples we have is more than 
          #         `effectively_independent_samples` (default 100) times.
          # - Specified convergence:
          #       + Burn for `burn` (default 2000) steps
          #       + Sample for `sample` (default 2000) steps
  
          kwd = {
              "auto_convergence": False, # TODO CHANGE ME
              "walkers": 100,
              "burn": 2000,
              "sample": 2000,
              # The minimum_sample, n_tau_exp_as_burn_in, minimum_eis are only
              # used if auto_convergence is turned on.
              "minimum_sample": 2000,
              "maximum_sample": 100000,
              "n_tau_exp_as_burn_in": 3,
              "minimum_effective_independent_samples": 100,
              "check_convergence_frequency": 1000,
              "a": 2.0,
              "threads": 1
          }
  
          # Update from the model, then update from any keyword arguments given.
          kwd.update(self._configuration.get("infer", {}).copy())
          kwd.update(**kwargs)
  
          # Make some checks.
          if kwd["walkers"] % 2 > 0 or kwd["walkers"] < 2 * len(parameters):
              raise ValueError("the number of walkers must be an even number and "
                  "be at least twice the number of model parameters")
  
          check_keywords = ["threads", "a"]
          if kwd["auto_convergence"]:
              logger.info("Convergence will be estimated automatically.")
              check_keywords += ["minimum_sample", "check_convergence_frequency",
                  "minimum_effective_independent_samples", "n_tau_exp_as_burn_in",
                  "maximum_sample"]
          
          else:
              check_keywords += ["burn", "sample"]
              logger.warn("No convergence checks will be done!")
              logger.info("Burning for {0} steps and sampling for {1} with {2} "\
                  "walkers".format(kwd["burn"], kwd["sample"], kwd["walkers"]))
  
          for keyword in check_keywords:
              if kwd[keyword] < 1:
                  raise ValueError("keyword {} must be a positive value".format(
                      keyword))
  
          # Check for non-standard proposal scales.
          if kwd["a"] != 2.0:
              logger.warn("Using proposal scale of {0:.2f}".format(kwd["a"]))
  
          # If no initial proposal given, estimate the model parameters.
          if initial_proposal is None:
              initial_proposal = self.estimate(data)
  
          # Initial proposal could be:
          #   - an array (N_walkers, N_dimensions)
          #   - a dictionary containing key/value pairs for the dimensions
          if isinstance(initial_proposal, dict):
  
              wavelengths_required = []
              for channel, spectrum in zip(matched_channels, data):
                  if channel is None: continue
                  z = initial_proposal.get("z",
                      initial_proposal.get("z_{}".format(channel), 0))
                  wavelengths_required.append(
                      [spectrum.disp[0] * (1 - z), spectrum.disp[-1] * (1 - z)])
  
              closest_point = [initial_proposal[p] \
                  for p in self.grid_points.dtype.names]
              subset_bounds = self._initialise_approximator(
                  closest_point=closest_point, 
                  wavelengths_required=wavelengths_required, force=True, **kwargs)
  
              initial_proposal = self._initial_proposal_distribution(
                  parameters, initial_proposal, kwd["walkers"])
  
          elif isinstance(initial_proposal, np.ndarray):
              initial_proposal = np.atleast_2d(initial_proposal)
              if initial_proposal.shape != (kwd["walkers"], len(parameters)):
                  raise ValueError("initial proposal must be an array of shape "\
                      "(N_parameters, N_walkers) ({0}, {1})".format(kwd["walkers"],
                          len(parameters)))
  
          # Prepare the convolution functions.
          self._create_convolution_functions(matched_channels, data, parameters)
  
          # Create the sampler.
          logger.info("Creating sampler with {0} walkers and {1} threads".format(
              kwd["walkers"], kwd["threads"]))
          debug = kwargs.get("debug", False)
          sampler = emcee.EnsembleSampler(kwd["walkers"], len(parameters),
              inference.ln_probability, a=kwd["a"], threads=kwd["threads"],
              args=(parameters, self, data, debug),
              kwargs={"matched_channels": matched_channels})
  
          # Regardless of whether we automatically check for convergence or not,
          # we will still need to burn in for some minimum amount of time.
          if kwd["auto_convergence"]:
              # Sample for `minimum_sample` period.
              descr, iterations = "", kwd["minimum_sample"]
          else:
              # Sample for `burn` period
              descr, iterations = "burn-in", kwd["burn"]
  
          # Start sampling.
          t_init = time()
          acceptance_fractions = []
          progress_bar = kwargs.get("__show_progress_bar", True)
          sampler, init_acceptance_fractions, pos, lnprob, rstate, init_elapsed \
              = self._sample(sampler, initial_proposal, iterations, descr=descr,
                  parameters=parameters, __show_progress_bar=progress_bar)
          acceptance_fractions.append(init_acceptance_fractions)
  
          # If we don't have to check for convergence, it's easy:
          if not kwd["auto_convergence"]:
  
              # Save the chain and log probabilities before we reset the chain.
              burn, sample = kwd["burn"], kwd["sample"]
              converged = None # we don't know!
              burn_chains = sampler.chain
              burn_ln_probabilities = sampler.lnprobability
  
              # Reset the chain.
              logger.debug("Resetting chain...")
              sampler.reset()
  
              # Sample the posterior.
              sampler, prod_acceptance_fractions, pos, lnprob, rstate, t_elapsed \
                  = self._sample(sampler, pos, kwd["sample"], lnprob0=lnprob,
                      rstate0=rstate, descr="production", parameters=parameters,
                      __show_progress_bar=progress_bar)
  
              production_chains = sampler.chain
              production_ln_probabilities = sampler.lnprobability
              acceptance_fractions.append(prod_acceptance_fractions)
  
          else:
  
              # Start checking for convergence at a frequency
              # of check_convergence_frequency
              last_state = [pos, lnprob, rstate]
              converged, total_steps = False, 0 + iterations
              min_eis_required = kwd["minimum_effective_independent_samples"]
              while not converged and kwd["maximum_sample"] > total_steps:
                  
                  # Check for convergence.
                  # Estimate the exponential autocorrelation time.
                  try:
                      tau_exp, rho, rho_max_fit \
                          = utils.estimate_tau_exp(sampler.chain)
  
                  except:
                      logger.exception("Exception occurred when trying to "
                          "estimate the exponential autocorrelation time:")
  
                      logger.info("To recover, we are temporarily setting tau_exp"
                          " to {0}".format(total_steps))
                      tau_exp = total_steps
  
                  logger.info("Estimated tau_exp at {0} is {1:.0f}".format(
                      total_steps, tau_exp))
  
                  # Grab everything n_tau_exp_as_burn_in times that.
                  burn = int(np.ceil(tau_exp)) * kwd["n_tau_exp_as_burn_in"]
                  sample = sampler.chain.shape[1] - burn
  
                  if 1 > sample:
                      logger.info("Sampler has not converged because {0}x the "
                          "estimated exponential autocorrelation time of {1:.0f}"
                          " is step {2}, and we are only at step {3}".format(
                              kwd["n_tau_exp_as_burn_in"], tau_exp, burn,
                              total_steps))
          
                  else:
  
                      # Calculate the integrated autocorrelation time in the 
                      # remaining sample, for every parameter.
                      tau_int = utils.estimate_tau_int(sampler.chain[:, burn:])
  
                      # Calculate the effective number of independent samples in 
                      # each parameter.
                      num_effective = (kwd["walkers"] * sample)/(2*tau_int)
                      
                      logger.info("Effective number of independent samples in "
                          "each parameter:")
                      for parameter, n_eis in zip(parameters, num_effective):
                          logger.info("\t{0}: {1:.0f}".format(parameter, n_eis))
  
                      if num_effective.min() > min_eis_required:
                          # Converged.
                          converged = True
                          logger.info("Convergence achieved ({0:.0f} > {1:.0f})"\
                              .format(num_effective.min() > min_eis_required))
  
                          # Separate the samples into burn and production..
                          burn_chains = sampler.chain[:, :burn, :]
                          burn_ln_probabilities = sampler.lnprobability[:burn]
  
                          production_chains = sampler.chain[:, burn:, :]
                          production_ln_probabilities = sampler.lnprobability[burn:]
                          break
  
                      else:
                          # Nope.
                          logger.info("Sampler has not converged because it did "
                              "not meet the minimum number of effective "
                              "independent samples ({0:.0f})".format(kwd["n"]))
                  
                  # Keep sampling.
                  iterations = kwd["check_convergence_frequency"]
                  logger.info("Trying for another {0} steps".format(iterations))
  
                  pos, lnprob, rstate = last_state
                  sampler, af, pos, lnprob, rstate, t_elapsed = self._sample(
                      sampler, pos, iterations, lnprob0=lnprob, rstate0=rstate,
                      descr="", parameters=parameters,
                      __show_progress_bar=progress_bar)
  
                  total_steps += iterations
                  acceptance_fractions.append(af)
                  last_state.extend(pos, lnprob, rstate)
                  del last_state[:3]
  
              if not converged:
                  logger.warn("Maximum number of samples ({:.0f}) reached without"
                      "convergence!".format(kwd["maximum_sample"]))
  
          logger.info("Total time elapsed: {0} seconds".format(time() - t_init))
          
          if sampler.pool:
              sampler.pool.close()
              sampler.pool.join()
  
          # Stack burn and production information together.
          chains = np.hstack([burn_chains, production_chains])
          lnprobability = np.hstack([
              burn_ln_probabilities, production_ln_probabilities])
          acceptance_fractions = np.hstack(acceptance_fractions)
  
          chi_sq, dof, model_fluxes = self._chi_sq(dict(zip(parameters, 
              [np.percentile(chains[:, burn:, i], 50) 
                  for i in range(len(parameters))])), data)
  
          # Convert velocity scales.
          symbol, scale, units = self._preferred_redshift_scale
          labels = [] + parameters
          scales = np.ones(len(parameters))
          if symbol != "z":
              for i, parameter in enumerate(parameters):
                  if parameter == "z" or parameter.startswith("z_"):
                      chains[:, :, i] *= scale
                      scales[i] = scale
                      if "_" in parameter:
                          labels[i] = "_".join([symbol, parameter.split("_")[1:]])
                      else:
                          labels[i] = symbol
                      logger.debug("Scaled {0} (now {1}) to units of {2}".format(
                          parameter, labels[i], units))
  
          # Calculate MAP values and associated uncertainties.
          theta = OrderedDict()
          for i, label in enumerate(labels):
              l, c, u = np.percentile(chains[:, burn:, i], [16, 50, 84])
              theta[label] = (c, u-c, l-c)
  
          # Re-arrange the chains to be in the same order as the model parameters.
          indices = np.array([parameters.index(p) \
              for p in self.parameters if p in parameters])
          chains = chains[:, :, indices]
  
          # Remove the convolution functions.
          if not kwargs.get("__keep_convolution_functions", False):
              self._destroy_convolution_functions()
  
          if full_output:
              metadata = {
                  "burn": burn,
                  "walkers": kwd["walkers"],
                  "sample": sample,
                  "parameters": labels,
                  "scales": scales,
                  "chi_sq": chi_sq,
                  "dof": dof
              }
              return (theta, chains, lnprobability, acceptance_fractions, sampler,
                  metadata)
          return theta