Potential numerical inaccuracy in weighted variance

[fweber144]

The identity Var(X) = E(X^2) - E(X)^2 is prone to numerical cancellation (see, e.g., here). In a similar form, this identity is used in the predictive variances from predvar_gauss() (perhaps at other places, too):

projpred/R/extend_family.R

Lines 125 to 126 in 6cec7d6

	mu_mean <- mu %*% wsample
	mu_var <- mu^2 %*% wsample - mu_mean^2

Perhaps this doesn't have any larger impact in projpred, but I wanted to highlight it anyway. Alternatives to the current implementation could be those mentioned here.