<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def _get_quantiles(self, X, width, quantiles, modelmat=None, lp=None, prediction=False, xform=True, term=-1):
 """ estimate prediction intervals for LinearGAM Parameters X : array input data of shape (n_samples, m_features) width : float on (0, 1) quantiles : array-like of floats on (0, 1) instead of specifying the prediciton width, one can specify the quantiles. so width=.95 is equivalent to quantiles=[.025, .975] modelmat : array of shape or None, default: None lp : array or None, default: None prediction : bool, default: True. whether to compute prediction intervals (True) or confidence intervals (False) xform : bool, default: True, whether to apply the inverse link function and return values on the scale of the distribution mean (True), or to keep on the linear predictor scale (False) term : int, default: -1 Returns ------- intervals: np.array of shape (n_samples, 2 or len(quantiles)) Notes ----- when the scale parameter is known, then we can proceed with a large sample approximation to the distribution of the model coefficients where B_hat ~ Normal(B, cov) when the scale parameter is unknown, then we have to account for the distribution of the estimated scale parameter, which is Chi-squared. since we scale our estimate of B_hat by the sqrt of estimated scale, we get a t distribution: Normal / sqrt(Chi-squared) ~ t see Simon Wood section 1.3.2, 1.3.3, 1.5.5, 2.1.5 """

        if quantiles is not None:
            quantiles = np.atleast_1d(quantiles)
        else:
            alpha = (1 - width)/2.
            quantiles = [alpha, 1 - alpha]
        for quantile in quantiles:
            if (quantile >= 1) or (quantile <= 0):
                raise ValueError('quantiles must be in (0, 1), but found {}'\
                                 .format(quantiles))

        if modelmat is None:
            modelmat = self._modelmat(X, term=term)
        if lp is None:
            lp = self._linear_predictor(modelmat=modelmat, term=term)

        idxs = self.terms.get_coef_indices(term)
        cov = self.statistics_['cov'][idxs][:, idxs]

        var = (modelmat.dot(cov) * modelmat.A).sum(axis=1)
        if prediction:
            var += self.distribution.scale

        lines = []
        for quantile in quantiles:
            if self.distribution._known_scale:
                q = sp.stats.norm.ppf(quantile)
            else:
                q = sp.stats.t.ppf(quantile, df=self.statistics_['n_samples'] -
                                                self.statistics_['edof'])

            lines.append(lp + q * var**0.5)
        lines = np.vstack(lines).T

        if xform:
            lines = self.link.mu(lines, self.distribution)
        return lines