private static int haasAndStokes(int nRows,int sampleSize,ReaderColumnSelection sampleRowsReader){
  HashMap<DblArray,Integer> valsCount=getValCounts(sampleRowsReader);
  if (valsCount.size() == 0)   return 1;
  int[] freqCounts=getFreqCounts(valsCount);
  float q=((float)sampleSize) / nRows;
  float _1MinusQ=1 - q;
  float duj1Fraction=((float)sampleSize) / (sampleSize - _1MinusQ * freqCounts[0]);
  float duj1=duj1Fraction * valsCount.size();
  float gamma=0;
  for (int i=1; i <= freqCounts.length; i++) {
    gamma+=i * (i - 1) * freqCounts[i - 1];
  }
  gamma*=duj1 / sampleSize / sampleSize;
  gamma+=duj1 / nRows - 1;
  gamma=Math.max(gamma,0);
  int estimate;
  if (gamma < HAAS_AND_STOKES_ALPHA1) {
    estimate=(int)(duj1Fraction * (valsCount.size() - freqCounts[0] * _1MinusQ * Math.log(_1MinusQ)* gamma / q));
  }
 else   if (gamma < HAAS_AND_STOKES_ALPHA2) {
    int numRemovedClasses=0;
    float updatedNumRows=nRows;
    int updatedSampleSize=sampleSize;
    for (    Integer cnt : valsCount.values()) {
      if (cnt > HAAS_AND_STOKES_UJ2A_C) {
        numRemovedClasses++;
        freqCounts[cnt - 1]--;
        updatedSampleSize-=cnt;
        updatedNumRows-=((float)cnt) / q;
      }
    }
    if (updatedSampleSize == 0) {
      estimate=(int)(duj1Fraction * (valsCount.size() - freqCounts[0] * (_1MinusQ) * Math.log(_1MinusQ)* gamma / q));
    }
 else {
      float updatedQ=((float)updatedSampleSize) / updatedNumRows;
      int updatedSampleCardinality=valsCount.size() - numRemovedClasses;
      float updatedDuj1Fraction=((float)updatedSampleSize) / (updatedSampleSize - (1 - updatedQ) * freqCounts[0]);
      float updatedDuj1=updatedDuj1Fraction * updatedSampleCardinality;
      float updatedGamma=0;
      for (int i=1; i <= freqCounts.length; i++) {
        updatedGamma+=i * (i - 1) * freqCounts[i - 1];
      }
      updatedGamma*=updatedDuj1 / updatedSampleSize / updatedSampleSize;
      updatedGamma+=updatedDuj1 / updatedNumRows - 1;
      updatedGamma=Math.max(updatedGamma,0);
      estimate=(int)(updatedDuj1Fraction * (updatedSampleCardinality - freqCounts[0] * (1 - updatedQ) * Math.log(1 - updatedQ)* updatedGamma / updatedQ)) + numRemovedClasses;
    }
  }
 else {
    float fraq1Numer=0;
    float fraq1Denom=0;
    float fraq2Numer=0;
    float fraq2Denom=0;
    for (int i=1; i <= freqCounts.length; i++) {
      fraq1Numer+=i * q * q* Math.pow(1 - q * q,i - 1)* freqCounts[i - 1];
      fraq1Denom+=Math.pow(_1MinusQ,i) * (Math.pow(1 + q,i) - 1) * freqCounts[i - 1];
      fraq2Numer+=Math.pow(_1MinusQ,i) * freqCounts[i - 1];
      fraq2Denom+=i * q * Math.pow(_1MinusQ,i - 1)* freqCounts[i - 1];
    }
    estimate=(int)(valsCount.size() + freqCounts[0] * fraq1Numer / fraq1Denom * fraq2Numer * fraq2Numer / fraq2Denom / fraq2Denom);
  }
  return estimate < 1 ? 1 : estimate;
}

private static int haasAndStokes ( int nRows , int sampleSize , ReaderColumnSelection sampleRowsReader ) { HashMap < DblArray , Integer > valsCount = getValCounts ( sampleRowsReader ) ; if ( valsCount . size ( ) == 0 ) return 1 ; int [ ] freqCounts = getFreqCounts ( valsCount ) ; float q = ( ( float ) sampleSize ) / nRows ; float _1MinusQ = 1 - q ; float duj1Fraction = ( ( float ) sampleSize ) / ( sampleSize - _1MinusQ * freqCounts [ 0 ] ) ; float duj1 = duj1Fraction * valsCount . size ( ) ; float gamma = 0 ; for ( int i = 1 ; i <= freqCounts . length ; i ++ ) { gamma += i * ( i - 1 ) * freqCounts [ i - 1 ] ; } gamma *= duj1 / sampleSize / sampleSize ; gamma += duj1 / nRows - 1 ; gamma = Math . max ( gamma , 0 ) ; int estimate ; if ( gamma < HAAS_AND_STOKES_ALPHA1 ) { estimate = ( int ) ( duj1Fraction * ( valsCount . size ( ) - freqCounts [ 0 ] * _1MinusQ * Math . log ( _1MinusQ ) * gamma / q ) ) ; } else if ( gamma < HAAS_AND_STOKES_ALPHA2 ) { int numRemovedClasses = 0 ; float updatedNumRows = nRows ; int updatedSampleSize = sampleSize ; for ( Integer cnt : valsCount . values ( ) ) { if ( cnt > HAAS_AND_STOKES_UJ2A_C ) { numRemovedClasses ++ ; freqCounts [ cnt - 1 ] -- ; updatedSampleSize -= cnt ; updatedNumRows -= ( ( float ) cnt ) / q ; } } if ( updatedSampleSize == 0 ) { estimate = ( int ) ( duj1Fraction * ( valsCount . size ( ) - freqCounts [ 0 ] * ( _1MinusQ ) * Math . log ( _1MinusQ ) * gamma / q ) ) ; } else { float updatedQ = ( ( float ) updatedSampleSize ) / updatedNumRows ; int updatedSampleCardinality = valsCount . size ( ) - numRemovedClasses ; float updatedDuj1Fraction = ( ( float ) updatedSampleSize ) / ( updatedSampleSize - ( 1 - updatedQ ) * freqCounts [ 0 ] ) ; float updatedDuj1 = updatedDuj1Fraction * updatedSampleCardinality ; float updatedGamma = 0 ; for ( int i = 1 ; i <= freqCounts . length ; i ++ ) { updatedGamma += i * ( i - 1 ) * freqCounts [ i - 1 ] ; } updatedGamma *= updatedDuj1 / updatedSampleSize / updatedSampleSize ; updatedGamma += updatedDuj1 / updatedNumRows - 1 ; updatedGamma = Math . max ( updatedGamma , 0 ) ; estimate = ( int ) ( updatedDuj1Fraction * ( updatedSampleCardinality - freqCounts [ 0 ] * ( 1 - updatedQ ) * Math . log ( 1 - updatedQ ) * updatedGamma / updatedQ ) ) + numRemovedClasses ; } } else { float fraq1Numer = 0 ; float fraq1Denom = 0 ; float fraq2Numer = 0 ; float fraq2Denom = 0 ; for ( int i = 1 ; i <= freqCounts . length ; i ++ ) { fraq1Numer += i * q * q * Math . pow ( 1 - q * q , i - 1 ) * freqCounts [ i - 1 ] ; fraq1Denom += Math . pow ( _1MinusQ , i ) * ( Math . pow ( 1 + q , i ) - 1 ) * freqCounts [ i - 1 ] ; fraq2Numer += Math . pow ( _1MinusQ , i ) * freqCounts [ i - 1 ] ; fraq2Denom += i * q * Math . pow ( _1MinusQ , i - 1 ) * freqCounts [ i - 1 ] ; } estimate = ( int ) ( valsCount . size ( ) + freqCounts [ 0 ] * fraq1Numer / fraq1Denom * fraq2Numer * fraq2Numer / fraq2Denom / fraq2Denom ) ; } return estimate < 1 ? 1 : estimate ; }
haas , peter j . , and lynne stokes . " estimating the number of classes in a finite population . " journal of the american statistical association 93 . 444 ( 1998 ) : 1475 - 1487 . the hybrid estimator given by eq . 33 in section 6
