def fit(
            self,
            data,
            states,
            estimator=None,
            complete_samples_only=True,
            **kwargs):
        """
        Determine βs from data

        Parameters
        ----------
        data: pandas.DataFrame
            Dataframe containing samples from the conditional distribution, p(Y|X)
            estimator: 'MLE' or 'MAP'

        completely_samples_only: boolean (True or False)
            Are they downsampled or complete? Defaults to True

        """
        if estimator == 'MLE':
            mean, variance = self.maximum_likelihood_estimator(data, states)
        elif estimator == 'MAP':
            raise NotImplementedError(
                "fit method has not been implemented using Maximum A-Priori (MAP)")

        return mean, variance