clipped statistics broken, multiple bugs

[Paul Price]

There are several bugs in the clipped statistics function that eliminate its
functionality. There exist no workarounds to these.

1. The function does not iterate the specified number of times: once a single

iteration is completed, the function gets sample statistics (mean, standard
deviation), and then if these are not NAN, immediately jumps out of the
iteration loop:

if (! (isnan(stats->sampleMean)

isnan(stats->sampleStdev))) {

clippedMean = stats->sampleMean;
clippedStdev = stats->sampleStdev;
Exit loop
i = stats->clipIter;

}

2. Values in the vector which are to be clipped are actually not clipped (in

addition to the above problem 1), due to transposing "i" for "j" in the mask:

for (j = 0; j < myVector->n; j++) {

if (fabs(myVector->data.F32[j] - clippedMean) >

(stats->clipSigma * errors->data.F32[j])) {

tmpMask->data.U8[i] = 0xff;

}

}

3. The function calculating the clipped mean uses the same psStats object to

calculate the sample statistics during the iterations, potentially overwriting
some data that the user wants to keep.

4. The clipped statistics function is actually run twice, since it appears in

two "if" statements, both of which are executed if the function does not return -1:

if ((stats->options & PS_STAT_CLIPPED_MEAN)

(stats->options &

PS_STAT_CLIPPED_STDEV)) {

if (-1 > p_psVectorClippedStats(inF32, errorsF32, mask, maskVal, stats)) {

...

} else if (-2 == p_psVectorClippedStats(inF32, errorsF32, mask, maskVal,

stats)) {

...

}

}

[Paul Price]

5. The masking of clipped values only takes effect if the maskVal specified by

the user is 0xff, since this is how the values are masked, and this is fed into
the sample mean calculator using the maskVal specified by the user:

tmpMask->data.U8[j] = 0xff;

[...]

p_psVectorSampleMean(myVector, errors, tmpMask, maskVal, stats);

[gusciora@…]

I fixed all these bugs except for the last one, which uses 0xff for the maskVal
during the iterations to mask clipped elements. I'm not convinced that code is
incorrect. The p_psVectorSampleMean() and p_psVectorSampleStdev() routines will
do a logical-and of the user-supplied maskVal with the tmpMask vector that we
are create. Therefore, any non-zero value for the user-supplied maskVal will
cause both the elements that the user wanted to be masked, and the elements that
we mask because they fall outside the clipped range, will be masked by those
functions. Let me know if this is incorrect.

Also, since we are on this topic, the ADD states that the sample median is used
as the first estimator of the mean, and then during the iterations, the sample
mean is used. This sounds plausible to me, but a little odd. Can you verify
that this is intended?

sampleMean and sampleStdev

[Paul Price]

As to bug 5, consider the case that the user has no mask of his own (the
clipping will do the masking for him). He simply wants the stats, so he calls:

stats = psStatsAlloc(PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
stats->clipSigma = 3.5;
stats->clipIter = 2;
(void)psVectorStats(stats, fluxes, NULL, NULL, 0);

Since he has no input mask, he figures it doesn't matter what maskVal he
specifies, and 0 is the default that immediately enters his head. His choice of
a value that he doesn't use shouldn't cause the function to behave incorrectly.

As to the question about the use of the median on the first iteration, that is
correct. The idea is that the median is less prone to outliers than the mean.
Once the clipping starts however, the mean should be good.

[gusciora@…]

You're right. I'll use 0xff for the maskVal when calling p_psVectorSampleMean()
and p_psVectorSampleStdev() inside that loop.

[Paul Price]

6. Prototype for psVectorStats is missing "const"s.

SDRS:
psStats *psVectorStats(psStats *stats,

const psVector *in,
const psVector *errors,
const psVector *mask,
unsigned int maskVal
);

psLib:
psStats* psVectorStats(

psStats* stats, /< stats structure defines stats to be calculated and how
psVector* in, /< Vector to be analysed.
psVector* errors, /< Errors.
psVector* mask, /< Ignore elements where (maskVector & maskVal) != 0:

must be INT or NULL

psU32 maskVal /< Only mask elements with one of these bits set in

maskVector
);

[Paul Price]

7. In the first iteration, the standard deviation (supposed to be the deviation

from the median for this iteration only) is calculated from the mean, not the
median. After adding some psTrace statements around the calls to
p_psVectorSampleMedian and p_psVectorSampleStdev, I get:

Start values: Mean: nan Median: nan StDev: nan
Pre-stdev values: Mean: nan Median: 52.142525 StDev: nan
Post-stdev values: Mean: -70.671822 Median: 52.142525 StDev: 1073.415039

8. tmpMask is not initialised to zeros if maskVector is NULL, meaning that (at

least on some architectures) there are values that are treated as masked even
though they haven't been explicitly masked.

[gusciora@…]

This should be done in the current version of the code (probably in the
April release as well, I think).

[Paul Price]

Bug has been resolved.... closing.

[Paul Price]

Bugs have been fixed... closing.