psStats robust stats not well designed

[eugene]

I've been looking at how we have spec-ed psStats robust statistics functions,
and I am not happy with some of the parameters. the resulting statistics are
somewhat sensitive to the choices of the bin size and smoothing scales (which
are also somewhat coupled) as well as the number of bins used in the fit (what
fraction of LQ to UQ is kept). I will explore this problem and consult with
people here to fix this algorithm.

Meanwhile, two improvements can be made:

• the clipped statistic is used to make a quick guess of MEAN & SIGMA. if the

dataset is large, this iteration can take some substantial amount of time. If
the number of points in the vector is larger than some number (let's set it at
50000, but use a #define), a random fraction of the pixels should be dropped to
keep the total number at about 50000 (this is only for the vector for which
clipped stats are measured).

• the fit of a gaussian to the points of interest in the peak should not use the

LM non-linear fitting. it is better to take the natural log of the ordinate and
fit a 2nd order polynomial (much faster). The resulting terms relate directly
to the gaussian parameters:

$SIGMA = sqrt (-0.5/$C2)
$MEAN = $C1 * $SIGMA^{2
$NORM = exp ($C0 + $MEAN}2 / (2*$SIGMA²⁾⁾

where C0, C1, C2 are the coefficients of the 0, 1, 2 order terms. the gaussian
which is fitted has the form:

f(x) = $NORM * exp(-0.5*(x - $MEAN)^{2 / $SIGMA}2)

[gusciora@…]

Okay. However, regarding the clipped stats, the latest ADD sets the cutoff for
data points as 10,000, while you request 50,000 in this bug. Status: not
implemented yet.

Regarding the fitting of the 2-D poly: not implemented yet. Currently, this
fails for other reasons (the change to the LMM minimization function prototype).

[gusciora@…]

Regarding the switch to a 2-D polynomial, how shall we determine the
coefficients? We could use psMinimizeLMChi2() with a 2-D polynomial instead of
a Gaussian. Or we could derive it in a similar manner to how we fit a
multi-dimensional polynomial to the sky background in the pilot project
(building and solving a matrix equation ...). Or is there some other method?

[gusciora@…]

I implemented the 1-D polynomial fit. The resulting numbers for the mean and
sigma do not agree with those obtained with the LMM fit. After I fit the
polynomial to the data, I calculate the sigma and mean based on your
instructions in this bug report:

psF32 polyFitSigma = PS_SQRT_F32(-0.5 / tmpPoly->coeff[2]);
psF32 polyFitMean = tmpPoly->coeff[1] * PS_SQR(polyFitSigma);

Can you review that?

[gusciora@…]

I implemented these

psF32 polyFitSigma = PS_SQRT_F32(-0.5 / tmpPoly->coeff[2]);
psF32 polyFitMean = tmpPoly->coeff[1] * PS_SQR(polyFitSigma);

but the resulting mean and sigma are not correct. Can you verify the
math is correct (possibly by looking at the code (approximately line 1644
of psStats.c in the correct CVS tree, maybe slightly different in the last
release of the code))?

[Paul Price]

* Bug 338 has been marked as a duplicate of this bug. *

[Paul Price]

I think this has been cleared up, and though the psStats code probably needs to
be cleaned some more, I think it does what it is supposed to.