public static double learnSymmetricConcentration(int[] countHistogram,
													 int[] observationLengths,
													 int numDimensions, 
													 double currentValue) {
		double currentDigamma;

		// The histogram arrays are presumably allocated before
		//  we knew what went in them. It is therefore likely that
		//  the largest non-zero value may be much closer to the 
		//  beginning than the end. We don't want to iterate over
		//  a whole bunch of zeros, so keep track of the last value.
		int largestNonZeroCount = 0;
		int[] nonZeroLengthIndex = new int[ observationLengths.length ];
		
		for (int index = 0; index < countHistogram.length; index++) {
			if (countHistogram[index] > 0) { largestNonZeroCount = index; }
		}

		int denseIndex = 0;
		for (int index = 0; index < observationLengths.length; index++) {
			if (observationLengths[index] > 0) {
				nonZeroLengthIndex[denseIndex] = index;
				denseIndex++;
			}
		}

		int denseIndexSize = denseIndex;

		for (int iteration = 1; iteration <= 200; iteration++) {
			
			double currentParameter = currentValue / numDimensions;

			// Calculate the numerator
			
			currentDigamma = 0;
			double numerator = 0;
		
			// Counts of 0 don't matter, so start with 1
			for (int index = 1; index <= largestNonZeroCount; index++) {
				currentDigamma += 1.0 / (currentParameter + index - 1);
				numerator += countHistogram[index] * currentDigamma;
			}
			
			// Now calculate the denominator, a sum over all observation lengths
			
			currentDigamma = 0;
			double denominator = 0;
			int previousLength = 0;
			
			double cachedDigamma = digamma(currentValue);
			
			for (denseIndex = 0; denseIndex < denseIndexSize; denseIndex++) {
				int length = nonZeroLengthIndex[denseIndex];
				
				if (length - previousLength > 20) {
					// If the next length is sufficiently far from the previous,
					//  it's faster to recalculate from scratch.
					currentDigamma = digamma(currentValue + length) - cachedDigamma;
				}
				else {
					// Otherwise iterate up. This looks slightly different
					//  from the previous version (no -1) because we're indexing differently.
					for (int index = previousLength; index < length; index++) {
						currentDigamma += 1.0 / (currentValue + index);
					}
				}
				
				denominator += currentDigamma * observationLengths[length];
			}
			
			currentValue = currentParameter * numerator / denominator;


			///System.out.println(currentValue + " = " + currentParameter + " * " + numerator + " / " + denominator);
		}

		return currentValue;
	}