def fit(self, min_ndata_factor=3, max_poly_order_factor=5, min_poly_order=2):
        """
        Fit the input data to the 'numerical eos', the equation of state employed
        in the quasiharmonic Debye model described in the paper:
        10.1103/PhysRevB.90.174107.

        credits: Cormac Toher

        Args:
            min_ndata_factor (int): parameter that controls the minimum number
                of data points that will be used for fitting.
                minimum number of data points =
                    total data points-2*min_ndata_factor
            max_poly_order_factor (int): parameter that limits the max order
                of the polynomial used for fitting.
                max_poly_order = number of data points used for fitting -
                                 max_poly_order_factor
            min_poly_order (int): minimum order of the polynomial to be
                considered for fitting.
        """
        warnings.simplefilter('ignore', np.RankWarning)

        get_rms = lambda x, y: np.sqrt(np.sum((np.array(x)-np.array(y))**2)/len(x))

        # list of (energy, volume) tuples
        e_v = [(i, j) for i, j in zip(self.energies, self.volumes)]
        ndata = len(e_v)
        # minimum number of data points used for fitting
        ndata_min = max(ndata - 2 * min_ndata_factor, min_poly_order + 1)
        rms_min = np.inf
        # number of data points available for fit in each iteration
        ndata_fit = ndata
        # store the fit polynomial coefficients and the rms in a dict,
        # where the key=(polynomial order, number of data points used for
        # fitting)
        all_coeffs = {}

        # sort by energy
        e_v = sorted(e_v, key=lambda x: x[0])
        # minimum energy tuple
        e_min = e_v[0]
        # sort by volume
        e_v = sorted(e_v, key=lambda x: x[1])
        # index of minimum energy tuple in the volume sorted list
        emin_idx = e_v.index(e_min)
        # the volume lower than the volume corresponding to minimum energy
        v_before = e_v[emin_idx - 1][1]
        # the volume higher than the volume corresponding to minimum energy
        v_after = e_v[emin_idx + 1][1]
        e_v_work = deepcopy(e_v)

        # loop over the data points.
        while (ndata_fit >= ndata_min) and (e_min in e_v_work):
            max_poly_order = ndata_fit - max_poly_order_factor
            e = [ei[0] for ei in e_v_work]
            v = [ei[1] for ei in e_v_work]
            # loop over polynomial order
            for i in range(min_poly_order, max_poly_order + 1):
                coeffs = np.polyfit(v, e, i)
                pder = np.polyder(coeffs)
                a = np.poly1d(pder)(v_before)
                b = np.poly1d(pder)(v_after)
                if a * b < 0:
                    rms = get_rms(e, np.poly1d(coeffs)(v))
                    rms_min = min(rms_min, rms * i / ndata_fit)
                    all_coeffs[(i, ndata_fit)] = [coeffs.tolist(), rms]
                    # store the fit coefficients small to large,
                    # i.e a0, a1, .. an
                    all_coeffs[(i, ndata_fit)][0].reverse()
            # remove 1 data point from each end.
            e_v_work.pop()
            e_v_work.pop(0)
            ndata_fit = len(e_v_work)

        logger.info("total number of polynomials: {}".format(len(all_coeffs)))

        norm = 0.
        fit_poly_order = ndata
        # weight average polynomial coefficients.
        weighted_avg_coeffs = np.zeros((fit_poly_order,))

        # combine all the filtered polynomial candidates to get the final fit.
        for k, v in all_coeffs.items():
            # weighted rms = rms * polynomial order / rms_min / ndata_fit
            weighted_rms = v[1] * k[0] / rms_min / k[1]
            weight = np.exp(-(weighted_rms ** 2))
            norm += weight
            coeffs = np.array(v[0])
            # pad the coefficient array with zeros
            coeffs = np.lib.pad(coeffs,
                                (0, max(fit_poly_order-len(coeffs), 0)),
                                'constant')
            weighted_avg_coeffs += weight * coeffs

        # normalization
        weighted_avg_coeffs /= norm
        weighted_avg_coeffs = weighted_avg_coeffs.tolist()
        # large to small(an, an-1, ..., a1, a0) as expected by np.poly1d
        weighted_avg_coeffs.reverse()

        self.eos_params = weighted_avg_coeffs
        self._set_params()