Maybe incorrect Brier score calculation

From https://github.com/tensorflow/tensorflow/issues/41664

System information

• Have I written custom code (as opposed to using a stock example script provided in TensorFlow): Yes
• TensorFlow version (use command below): 1.15.0
• TensorFlow Probability version: 0.8.0

Describe the current behavior

I believe that the way the Brier score calculation is implemented in tensorflow_probability.stats.brier_score is incorrect. The formula used is

sum_i p[i]* p[i] - 2*p[k]

where p is the probability vector over all discrete outcomes, and k is the realized outcome. (Note: This gives element wise Brier scores, and to get the actual Brier score across the dataset, reduce_mean() must be called.)

Describe the expected behavior

This formula does not match the reference cited nor Wikipedia nor the definition that sklearn uses. Also, it is stated in the tensorflow_probability docs that the Brier score can be negative, which is not true.

The reference cited is Brier’s original paper, found here, which states the formula as

where r is the number of possible classes, n is the number of forecasts, fij is the probability forecast of class j for instance i, and Eij is 0 or 1 depending if the event occurred in class j or not. This is the same formula that is used in sklearn and Wikipedia. By definition, this score cannot be negative.

If we consider the element wise Brier score from this formu