def deviance_residuals(self, X, y, weights=None, scaled=False):
        """
        method to compute the deviance residuals of the model

        these are analogous to the residuals of an OLS.

        Parameters
        ----------
        X : array-like
            Input data array of shape (n_saples, m_features)
        y : array-like
            Output data vector of shape (n_samples,)
        weights : array-like shape (n_samples,) or None, optional
            Sample weights.
            if None, defaults to array of ones
        scaled : bool, optional
            whether to scale the deviance by the (estimated) distribution scale

        Returns
        -------
        deviance_residuals : np.array
            with shape (n_samples,)
        """
        if not self._is_fitted:
            raise AttributeError('GAM has not been fitted. Call fit first.')

        y = check_y(y, self.link, self.distribution, verbose=self.verbose)
        X = check_X(X, n_feats=self.statistics_['m_features'],
                    edge_knots=self.edge_knots_, dtypes=self.dtype,
                    features=self.feature, verbose=self.verbose)
        check_X_y(X, y)

        if weights is not None:
            weights = np.array(weights).astype('f').ravel()
            weights = check_array(weights, name='sample weights',
                                  ndim=1, verbose=self.verbose)
            check_lengths(y, weights)
        else:
            weights = np.ones_like(y).astype('float64')

        mu = self.predict_mu(X)
        sign = np.sign(y-mu)
        return sign * self.distribution.deviance(y, mu,
                                                 weights=weights,
                                                 scaled=scaled) ** 0.5