obs_diag rmse normalization

[timhoar]

obs_diag computes the root-mean-squared-error (rmse) which begs the question: 'Error from what?'
The observation is considered to be the truth. This has subtle ramifications when normalizing the rmse ... do you normalize by 'N' or 'N-1' ? The reason statisticians normalize by N-1 is to account for the loss of the degree of freedom when you have to use the data to estimate the 'truth' used to compute the error. Since obs_diag presumes the observation is the truth, it has been normalizing by N.

[kdraeder]

I'm not a statistician, but normalizing by N sounds justified.
In addition, normalizing by N-1 causes uncomfortable downstream consequences.
If there's only 1 observation in a region at time, but that happens repeatedly
in some time span, then the RMSE of each of those is assigned NaN,
since division by 0 is not defined. Those NaNs are written to obs_diag_output.nc.
Then the plot_rmse_XXX_evolution script reads them in and tries to calculate
the grand ensemble statistics. Those are calculated over all of the times,
so there are > 1 observations in the grand ensemble. But they are all NaNs,
so the grand ensemble RMSE is printed on the figures as NaN,
even though there are > 1 observations, which were all used.
Our view is that people shouldn't be looking at plots with at most 1 ob per time,
but people probaby will, and may pester us about this apparent inconsistency.
If we're happy to use this "teaching moment" whenever it comes up,
then normalizing by N-1 is fine. Otherwise, N is the better choice.

[kdraeder]

Another implication of excluding single observation times from the grand ensemble statistics
is that in cases where there are times with > 1 ob, and numbers (instead of NaNs)
are reported for the statistics, those numbers will not be correct because the single ob
times were excluded. This error will be unimportant for large grand ensembles,
but there's no way for the average viewer to know how big the effect is for a given ensemble size.
Even if they know to count the ignored times and compare that with the whole ensemble size,
they can't know the (possibly nonlinear) effect that including those observations would have.
It's likely that isolated obs have larger biases and RMSEs than the ensemble as a whole
precisely because they are isolated; the model state may not have been pushed much to match
the observations.