variance/sd over 0-vector or 1-vector?

[Bob Carpenter]

We're using (n-1) in the denominator for
variance and standard deviation, but I can't
recall what we decided the boundary behavior
should be.

As the code stands now, we throw a domain
error in variance() and sd() for any input
with n < 2.

I recall discussing that behavior, but I never
recorded our decision. For n=0, the limiting
behavior is 0/-1, or 0. For n=1, the limiting
behavior is 0/0, or ???.

I expressed concern about trying to return a
finite value for n=0 or n=1 because it'll propagate
through other uses without raising an error
anywhere.

(The context of the question is that I'm adding bounds
checks to all the matrix ops for 1.0.2 to avoid seg
faulting in Eigen and our own libs. In so doing, I'm
looking closely again at matrix/vector op error checking.)

- Bob

[Andrew Gelman]

I think we should report 0 for variance and sd when n=0 or n=1. In any reasonable weighting scheme, variance estimates based on 0 or 1 point would be zero-weighted anyway, but I don't see a reason for crashing the program or sending out NA's just because we happen to compute a variance of a vector of length 1.
A

> --
>

[Bob Carpenter]

On 9/21/12 6:07 PM, Andrew Gelman wrote:
> I think we should report 0 for variance and sd when n=0 or n=1.
> In any reasonable weighting scheme, variance estimates based on 0 or 1
> point would be zero-weighted anyway,

This is all coming back to me, and I'm wondering
if we're talking about the same situation.

What do you mean by "report"? I'm not talking
about ESS or R-hat here, but the functions that
get called in Stan models.

And what's a "weighting scheme"? I'm used to
seeing weighting by precision (inverse variance).

I'm worried about catching cases where programmers
make an error and don't know it.

My limited imagination can't dream up a situation where
a user would want a return of 0 for size-0 or size-1
inputs.

I can just do what I'm told without understaning,
so if you really think it's important, I'll put it
in and document it without understanding why we
do it the way you're suggesting.

The main use I've seen so far is data normalization.

vector[K] y_raw;
...
vector[K] y;
...
real y_sd;

y_sd <- sd(y_raw);
y_mean <- mean(y_raw);
y <- (y_raw - y_mean) / y_sd; // I think we vectorized - scalar

In this case, if K = 1, then setting sd(y_raw)
to 0 sets the value in y[0] to 0 / 0, which
yields NaN.

In general, it's better to fail early as soon
as the problem comes up (here computing sd()),
because it's easier to debug.

This also brings up what to do with min, max, and
mean in the face of zero-length inputs. We currently
are raising exceptions if any of these is fed
a length-zero input.

> but I don't see a reason for crashing
> the program or sending out NA's
> just because we happen to compute a variance of a vector of length 1.

Nobody is suggesting either crashing or returning
NA. At most, we'd raise a domain error in the
same way we would if someone tries to supply a negative
scale parameter to the normal distribution.

NA values are an R thing for representing undefined
values. There's no NA in C++. We do use
NaN values depending on the error policy
configuration.




- Bob

[Andrew Gelman]

Bob,
I have no experience with zero-length arrays so I have no idea what these functions should return in that case. But for arrays of length 1, yes, I think that setting variance=0 is reasonable.
Consider the following calculation:
Sum_j (n_j-1)*var(y_j),
where each "j" is a vector y_j of length n_j.
In that case, we want "(n_j-1)*var(y_j)" to be 0 when n_j=1. Setting var(y_j) to zero for vectors of length 1 will do this.
A

> --
>

[Bob Carpenter]

OK, I'll set sd() and var() up to only throw an exceptions
for length 0 arrays and to return 1 for size 1 arrays and vectors.

If x = { x1 } is a size 1 array, then

var(x) = (x1 - mean(x))^2 / (1 - size(x))
= (x1 - x1)^2 / (1 - 1)
= 0 / 0

My thinking for making it throw a domain error is that
the IEEE arithmetic behavior for 0 / 0 is to return NaN.
But we can make it 1, as I said above.

Unless anyone objects, I'll also go with the usual
mathematical limits for min and max:

min(y) = +inf
max(y) = -inf

That is, min of an empty list returns the biggest
value and max returns the smallest value. The motivation
is that it does the right thing under array concatenation:

min(min(y),min(z)) = min(y+z)

so that if y is empty so that y+z = z, then

min(min({}),min(z)) = min(+inf,min(z)) = min(z)

- Bob

P.S. Could someone please just humor this poor
computer scientist and tell him why you want
(n - 1) * var(y) = 0 when y is an array of
length n = 1. Any finite value for var(y)
works because (n - 1) = 0. ???

[Jiqiang Guo]

Is the following 1 a typo?  Bot 0 and NaN are reasonable to me.  

Jiqiang

[Bob Carpenter]

On 9/22/12 3:48 PM, Jiqiang Guo wrote:
> Is the following 1 a typo? Bot 0 and NaN are reasonable to me.
>
> Jiqiang
>

> On Sat, Sep 22, 2012 at 2:11 AM, Bob Carpenter <ca...@alias-i.com <mailto:ca...@alias-i.com>> wrote:
>
> OK, I'll set sd() and var() up to only throw an exceptions
> for length 0 arrays and to return 1 for size 1 arrays and vectors.

^^ typo!

Yes, I meant 0. I still don't understand why, but it
really doesn't matter other than for the doc.

- Bob