docs: def maximum_likelihood_estimator(self, data, states):
          """ 
          Fit using MLE method.
  
          Parameters
          ----------
          data: pandas.DataFrame or 2D array
              Dataframe of values containing samples from the conditional distribution, (Y|X)
              and corresponding X values.
  
          states: All the input states that are jointly gaussian.
  
          Returns
          -------
          beta, variance (tuple): Returns estimated betas and the variance.
          """
          x_df = pd.DataFrame(data, columns=states)
          x_len = len(self.evidence)
  
          sym_coefs = []
          for i in range(0, x_len):
              sym_coefs.append('b' + str(i + 1) + '_coef')
  
          sum_x = x_df.sum()
          x = [sum_x['(Y|X)']]
          coef_matrix = pd.DataFrame(columns=sym_coefs)
  
          # First we compute just the coefficients of beta_1 to beta_N.
          # Later we compute beta_0 and append it.
          for i in range(0, x_len):
              x.append(self.sum_of_product(x_df['(Y|X)'], x_df[self.evidence[i]]))
              for j in range(0, x_len):
                  coef_matrix.loc[i, sym_coefs[j]] = self.sum_of_product(
                      x_df[self.evidence[i]], x_df[self.evidence[j]])
  
          coef_matrix.insert(0, 'b0_coef', sum_x[self.evidence].values)
          row_1 = np.append([len(x_df)], sum_x[self.evidence].values)
          coef_matrix.loc[-1] = row_1
          coef_matrix.index = coef_matrix.index + 1  # shifting index
          coef_matrix.sort_index(inplace=True)
  
          beta_coef_matrix = np.matrix(coef_matrix.values, dtype='float')
          coef_inv = np.linalg.inv(beta_coef_matrix)
          beta_est = np.array(np.matmul(coef_inv, np.transpose(x)))
          self.beta = beta_est[0]
  
          sigma_est = 0
          x_len_df = len(x_df)
          for i in range(0, x_len):
              for j in range(0, x_len):
                  sigma_est += self.beta[i + 1] * self.beta[j + 1] * (self.sum_of_product(
                      x_df[self.evidence[i]], x_df[self.evidence[j]]) / x_len_df - np.mean(x_df[self.evidence[i]]) * np.mean(x_df[self.evidence[j]]))
  
          sigma_est = np.sqrt(
              self.sum_of_product(
                  x_df['(Y|X)'],
                  x_df['(Y|X)']) /
              x_len_df -
              np.mean(
                  x_df['(Y|X)']) *
              np.mean(
                  x_df['(Y|X)']) -
              sigma_est)
          self.sigma_yx = sigma_est
          return self.beta, self.sigma_yx