private double[] findAndEvaluateThreshold(DoubleDynamicHistogram histogram){
  int n=histogram.getNumBins();
  double[] p1=new double[n];
  double[] p2=new double[n];
  double[] mu1=new double[n];
  double[] mu2=new double[n];
  double[] sigma1=new double[n];
  double[] sigma2=new double[n];
  double[] jt=new double[n];
{
    MeanVariance mv=new MeanVariance();
    DoubleHistogram.Iter forward=histogram.iter();
    for (int i=0; forward.valid(); i++, forward.advance()) {
      p1[i]=forward.getValue() + ((i > 0) ? p1[i - 1] : 0);
      mv.put(i,forward.getValue());
      mu1[i]=mv.getMean();
      sigma1[i]=mv.getNaiveStddev();
    }
  }
{
    MeanVariance mv=new MeanVariance();
    DoubleHistogram.Iter backwards=histogram.iter();
    backwards.seek(histogram.getNumBins() - 1);
    for (int j=n - 1; backwards.valid(); j--, backwards.retract()) {
      p2[j]=backwards.getValue() + ((j + 1 < n) ? p2[j + 1] : 0);
      mv.put(j,backwards.getValue());
      mu2[j]=mv.getMean();
      sigma2[j]=mv.getNaiveStddev();
    }
  }
  for (int i=0; i < n; i++) {
    jt[i]=1.0 + 2 * (p1[i] * (Math.log(sigma1[i]) - Math.log(p1[i])) + p2[i] * (Math.log(sigma2[i]) - Math.log(p2[i])));
  }
  int bestpos=-1;
  double bestgoodness=Double.NEGATIVE_INFINITY;
  double devPrev=jt[1] - jt[0];
  for (int i=1; i < jt.length - 1; i++) {
    double devCur=jt[i + 1] - jt[i];
    if (devCur >= 0 && devPrev <= 0) {
      double lowestMaxima=Double.POSITIVE_INFINITY;
      for (int j=i - 1; j > 0; j--) {
        if (jt[j - 1] < jt[j]) {
          lowestMaxima=Math.min(lowestMaxima,jt[j]);
          break;
        }
      }
      for (int j=i + 1; j < n - 2; j++) {
        if (jt[j + 1] < jt[j]) {
          lowestMaxima=Math.min(lowestMaxima,jt[j]);
          break;
        }
      }
      double localDepth=lowestMaxima - jt[i];
      final double mud=mu1[i] - mu2[i];
      double discriminability=mud * mud / (sigma1[i] * sigma1[i] + sigma2[i] * sigma2[i]);
      if (Double.isNaN(discriminability)) {
        discriminability=-1;
      }
      double goodness=localDepth * discriminability;
      if (goodness > bestgoodness) {
        bestgoodness=goodness;
        bestpos=i;
      }
    }
    devPrev=devCur;
  }
  Iter iter=histogram.iter();
  iter.seek(bestpos);
  return new double[]{iter.getRight(),bestgoodness};
}

private double [ ] findAndEvaluateThreshold ( DoubleDynamicHistogram histogram ) { int n = histogram . getNumBins ( ) ; double [ ] p1 = new double [ n ] ; double [ ] p2 = new double [ n ] ; double [ ] mu1 = new double [ n ] ; double [ ] mu2 = new double [ n ] ; double [ ] sigma1 = new double [ n ] ; double [ ] sigma2 = new double [ n ] ; double [ ] jt = new double [ n ] ; { MeanVariance mv = new MeanVariance ( ) ; DoubleHistogram . Iter forward = histogram . iter ( ) ; for ( int i = 0 ; forward . valid ( ) ; i ++ , forward . advance ( ) ) { p1 [ i ] = forward . getValue ( ) + ( ( i > 0 ) ? p1 [ i - 1 ] : 0 ) ; mv . put ( i , forward . getValue ( ) ) ; mu1 [ i ] = mv . getMean ( ) ; sigma1 [ i ] = mv . getNaiveStddev ( ) ; } } { MeanVariance mv = new MeanVariance ( ) ; DoubleHistogram . Iter backwards = histogram . iter ( ) ; backwards . seek ( histogram . getNumBins ( ) - 1 ) ; for ( int j = n - 1 ; backwards . valid ( ) ; j -- , backwards . retract ( ) ) { p2 [ j ] = backwards . getValue ( ) + ( ( j + 1 < n ) ? p2 [ j + 1 ] : 0 ) ; mv . put ( j , backwards . getValue ( ) ) ; mu2 [ j ] = mv . getMean ( ) ; sigma2 [ j ] = mv . getNaiveStddev ( ) ; } } for ( int i = 0 ; i < n ; i ++ ) { jt [ i ] = 1.0 + 2 * ( p1 [ i ] * ( Math . log ( sigma1 [ i ] ) - Math . log ( p1 [ i ] ) ) + p2 [ i ] * ( Math . log ( sigma2 [ i ] ) - Math . log ( p2 [ i ] ) ) ) ; } int bestpos = - 1 ; double bestgoodness = Double . NEGATIVE_INFINITY ; double devPrev = jt [ 1 ] - jt [ 0 ] ; for ( int i = 1 ; i < jt . length - 1 ; i ++ ) { double devCur = jt [ i + 1 ] - jt [ i ] ; if ( devCur >= 0 && devPrev <= 0 ) { double lowestMaxima = Double . POSITIVE_INFINITY ; for ( int j = i - 1 ; j > 0 ; j -- ) { if ( jt [ j - 1 ] < jt [ j ] ) { lowestMaxima = Math . min ( lowestMaxima , jt [ j ] ) ; break ; } } for ( int j = i + 1 ; j < n - 2 ; j ++ ) { if ( jt [ j + 1 ] < jt [ j ] ) { lowestMaxima = Math . min ( lowestMaxima , jt [ j ] ) ; break ; } } double localDepth = lowestMaxima - jt [ i ] ; final double mud = mu1 [ i ] - mu2 [ i ] ; double discriminability = mud * mud / ( sigma1 [ i ] * sigma1 [ i ] + sigma2 [ i ] * sigma2 [ i ] ) ; if ( Double . isNaN ( discriminability ) ) { discriminability = - 1 ; } double goodness = localDepth * discriminability ; if ( goodness > bestgoodness ) { bestgoodness = goodness ; bestpos = i ; } } devPrev = devCur ; } Iter iter = histogram . iter ( ) ; iter . seek ( bestpos ) ; return new double [ ] { iter . getRight ( ) , bestgoodness } ; }
evaluate the histogram to find a suitable threshold
