<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def _estimate_model_statistics(self, y, modelmat, inner=None, BW=None, B=None, weights=None, U1=None):
 """ method to compute all of the model statistics results are stored in the 'statistics_' attribute of the model, as a dictionary keyed by: - edof: estimated degrees freedom - scale: distribution scale, if applicable - cov: coefficient covariances - se: standarrd errors - AIC: Akaike Information Criterion - AICc: corrected Akaike Information Criterion - pseudo_r2: dict of Pseudo R-squared metrics - GCV: generailized cross-validation or - UBRE: Un-Biased Risk Estimator - n_samples: number of samples used in estimation Parameters y : array-like output data vector of shape (n_samples,) modelmat : array-like, default: None contains the spline basis for each feature evaluated at the input inner : array of intermediate computations from naive optimization BW : array of intermediate computations from either optimization B : array of intermediate computations from stable optimization weights : array-like shape (n_samples,) or None, default: None containing sample weights U1 : cropped U matrix from SVD. Returns ------- None """

        lp = self._linear_predictor(modelmat=modelmat)
        mu = self.link.mu(lp, self.distribution)
        self.statistics_['edof_per_coef'] = np.diagonal(U1.dot(U1.T))
        self.statistics_['edof'] = self.statistics_['edof_per_coef'].sum()
        if not self.distribution._known_scale:
            self.distribution.scale = self.distribution.phi(y=y, mu=mu, edof=self.statistics_['edof'], weights=weights)
        self.statistics_['scale'] = self.distribution.scale
        self.statistics_['cov'] = (B.dot(B.T)) * self.distribution.scale # parameter covariances. no need to remove a W because we are using W^2. Wood pg 184
        self.statistics_['se'] = self.statistics_['cov'].diagonal()**0.5
        self.statistics_['AIC'] = self._estimate_AIC(y=y, mu=mu, weights=weights)
        self.statistics_['AICc'] = self._estimate_AICc(y=y, mu=mu, weights=weights)
        self.statistics_['pseudo_r2'] = self._estimate_r2(y=y, mu=mu, weights=weights)
        self.statistics_['GCV'], self.statistics_['UBRE'] = self._estimate_GCV_UBRE(modelmat=modelmat, y=y, weights=weights)
        self.statistics_['loglikelihood'] = self._loglikelihood(y, mu, weights=weights)
        self.statistics_['deviance'] = self.distribution.deviance(y=y, mu=mu, weights=weights).sum()
        self.statistics_['p_values'] = self._estimate_p_values()