6
$\begingroup$

I have a data set which fits a logNormal distribution quite well. (From a theoretical point of view, it is some hard-to-tackle quotient distribution).

However, the data is quite dirty, so parameter estimation is far from trivial.

Right now, my approach is the following:

  1. Shift the distribution such that the minimum is almost 0.
  2. Logspace the data
  3. Use a robust Median and MAD parameter estimation (see Estimating parameters of a normal distribution: median instead of mean? for details)

The result is significantly better than before (Maximum difference from empirical CDF 0.034 instead of 0.081 and 0.224 without using MAD). It's not perfect in particular on the long tail where I expect outliers. The additional location parameter helped a lot. However, using the minimum is a very crude heuristic. I obviously cannot expect to observe the true minimum, but depending on the sample size the observed minimum will always be some small x larger.

Do you know any robust parameter estimation method (+ a reference if possible) for the $e^{\mathcal{N}(\mu, \sigma)} + c$ distribution family?

Note that e.g. scipy.stats.lognorm does also have such an additional third location parameter, just like the one I'm using, but I'm working in Java with my own code.

Update: I've just come across a thesis on this topic:

  • Estimating the Parameters of the Three-Parameter Lognormal Distribution
    Rodrigo J. Aristizabal

which includes pointers to some relevant literature, in particular to

  • Estimating Parameters of Logarithmic-Normal Distributions by Maximum Likelihood
    A. C. Cohen, Jr.

but I find it hard to get a formula out of these publications that I could implement.

$\endgroup$
5
  • $\begingroup$ @user603 I'm actually working in outlier detection. I need the robust estimations for advanced outlier detection outside of the "larger than 3 sigma" ideas, but e.g. on ensemble techniques for outlier detection. $\endgroup$ – Erich Schubert Feb 7 '13 at 16:34
  • $\begingroup$ @user603 I have also tried Levenberg Marquard fitting to the KDE curve of a truncated part of the data set, but that does not perform as good as the robust statistics directly. $\endgroup$ – Erich Schubert Feb 7 '13 at 16:35
  • $\begingroup$ In particular, I'm not interested in a binary "yes/no" outlier decision, but into judging and calibrating the "outlierness" of objects. $\endgroup$ – Erich Schubert Feb 7 '13 at 16:46
  • $\begingroup$ possible duplicate of Robust parameter estimation for Exponentially modified Gaussian distribution $\endgroup$ – Ellie Kesselman Aug 15 '14 at 5:32
  • $\begingroup$ @FeralOink that question is about a quite different distribution. Just as the uniform distribution is not the same as the normal distribution. $\endgroup$ – Erich Schubert Aug 18 '14 at 7:55
2
$\begingroup$

In case anyone is still interested, I have managed to implement Aristizabal's formulae in Java. This is more proof-of-concept than the requested "robust" code, but it is a starting point.

/**
 * Computes the point estimate of the shift offset (gamma) from the given sample. The sample array will be sorted by this method.<p>
 * Cf. Aristizabal section 2.2 ff.
 * @param sample {@code double[]}, will be sorted
 * @return gamma point estimate
 */
public static double pointEstimateOfGammaFromSample(double[] sample) {
    Arrays.sort(sample);
    DoubleUnaryOperator func = x->calculatePivotalOfSortedSample(sample, x)-1.0;
    double upperLimit = sample[0];
    double lowerLimit = 0;
    double gamma = bisect(func, lowerLimit, upperLimit);
    return gamma;
}

/**
 * Cf. Aristizabal's equation (2.3.1)
 * @param sample {@code double[]}, should be sorted in ascending order
 * @param gamma shift offset
 * @return pivotal value of sample
 */
private static double calculatePivotalOfSortedSample(final double[] sample, double gamma) {
    final int n=sample.length;
    final int n3=n/3;
    final double mid = avg(sample, gamma, n3+1, n-n3);
    final double low = avg(sample, gamma, 1, n3);
    final double upp = avg(sample, gamma, n-n3+1, n);
    final double result = (mid-low)/(upp-mid);
    return result;
}

/**
 * Computes average of sample values from {@code sample[l-1]} to {@code sample[u-1]}.
 * @param sample {@code double[]}, should be sorted in ascending order
 * @param gamma shift offset
 * @param l lower limit
 * @param u upper limit
 * @return average
 */
private static double avg(double[] sample, double gamma, int l, int u) {
    double sum = 0.0;
    for (int i=l-1;i<u;sum+=Math.log(sample[i++]-gamma));
    final int n = u-l+1;
    return sum/n;
}

/**
 * Naive bisection implementation. Should always complete if the given values actually straddles the root.
 * Will call {@link #secant(DoubleUnaryOperator, double, double)} if they do not, in which case the
 * call may not complete.
 * @param func Function solve for root value
 * @param lowerLimit Some value for which the given function evaluates < 0
 * @param upperLimit Some value for which the given function evaluates > 0
 * @return x value, somewhere between the lower and upper limits, which evaluates close enough to zero
 */
private static double bisect(DoubleUnaryOperator func, double lowerLimit, double upperLimit) {
    final double eps = 0.000001;
    double low=lowerLimit;
    double valAtLow = func.applyAsDouble(low);
    double upp=upperLimit;
    double valAtUpp = func.applyAsDouble(upp);
    if (valAtLow*valAtLow>0) {
        // Switch to secant method
        return secant(func, lowerLimit, upperLimit);
    }
    System.out.printf("bisect %f@%f -- %f@%f%n", valAtLow, low, valAtUpp, upp);
    double mid;
    while(true) {
        mid = (upp+low)/2;
        if (Math.abs(upp-low)/low<eps)
            break;
        double val = func.applyAsDouble(mid);
        if (Math.abs(val)<eps)
            break;
        if (val<0)
            low=mid;
        else
            upp=mid;
    }
    return mid;
}

/**
 * Naive secant root solver implementation. May not complete if root not found.
 * @param f Function solve for root value
 * @param a Some value for which the given function evaluates
 * @param b Some value for which the given function evaluates
 * @return x value which evaluates close enough to zero
 */
static double secant(final DoubleUnaryOperator f, double a, double b) {
    double fa = f.applyAsDouble(a);
    if (fa==0)
        return a;
    double fb = f.applyAsDouble(b);
    if (fb==0)
        return b;
    System.out.printf("secant %f@%f -- %f@%f%n", fa, a, fb, b);
    if (fa*fb<0) {
        return bisect(f, a, b);
    }
    while ( abs(b-a) > abs(0.00001*a) ) {
          final double m = (a+b)/2;
          final double k = (fb-fa)/(b-a);
          final double fm = f.applyAsDouble(m);
          final double x = m-fm/k;
          if (Math.abs(fa)<Math.abs(fb)) {
              // f(a)<f(b); Choose x and a
              b=x;
              fb=f.applyAsDouble(b);
          } else {
              // f(a)>=f(b); Choose x and b
              a=x;
              fa=f.applyAsDouble(a);
          }
          if (fa==0)
              return a;
          if (fb==0)
              return b;
          if (fa*fb<0) {
              // Straddling root; switch to bisect method
              return bisect(f, a, b);
          }
      }
    return (a+b)/2;

}
$\endgroup$
4
  • 1
    $\begingroup$ ELKI has L-Moment based estimators now, these are rather robust: elki.dbs.ifi.lmu.de/releases/release0.7.1/doc/de/lmu/ifi/dbs/… and IIRC can also estimate location. $\endgroup$ – Erich Schubert Sep 8 '16 at 7:55
  • $\begingroup$ Thanks, lenborge. For some reason I can't seem to be able to run your code. Here is my example: public static void main(String[] args) { double[] sample = {5.883585, 8.083535, 6.712414, 5.654602, 9.265344, 8.666495, 6.322371, 7.058905, 5.836143, 6.627015, 7.673438, 8.989364, 8.174420, 8.207641, 5.630372, 7.315935, 6.441540, 8.277943, 12.792464, 8.216356}; System.out.println(pointEstimateOfGammaFromSample(sample)); }` The program outputs bisect 0.042903@0.000000 -- Infinity@5.630372 and then goes in an infinite loop. $\endgroup$ – dragice Jan 10 '17 at 13:45
  • $\begingroup$ The "bisect" method, which is used to find a root of a function, requires that the root lies somewhere within the range of values given. You could either use some other numerical solver (I just used the bisect method as a quick-and-dirty way forward), or add some dummy data point with a large enough value. I tried just adding "100.0" to your data sample, and then the algorithm completed. $\endgroup$ – lenborje Jan 11 '17 at 14:42
  • $\begingroup$ If edited the code and added a secant solver which is used if the values given to bisect do not straddle the root. The variant should run with your sample data. $\endgroup$ – lenborje Jan 12 '17 at 9:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.