self.beta_estimate = multi_dot([self.inverse_covariance, X_ones.T, y])
// now we need the estimate of the noise variance
// reference: https://stat.ethz.ch/R-manual/R-devel/library/stats/html/summary.lm.html
residuals = (y - self.predict(X))
// square all the residuals
residuals **= 2
self.sigma_squared_estimate = residuals.sum() / max((n - d), 1)
self.covar = self.sigma_squared_estimate * self.inverse_covariance
After Change
// get the residual of the predictions and square it
pred -= y
pred **= 2
sum_squared_residuals = pred.sum()
self.sigma_squared_estimate = sum_squared_residuals / max((n - d), 1)
self.covar = self.sigma_squared_estimate * self.inverse_covariance
def predict(self, X, random_draw=False):