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
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,