Java Code Examples for org.apache.commons.math3.special.Erf#erfc()

/**
 * returns the probability of x falling within the range of a to b
 * under a normal distribution with mean mu and standard deviation sigma
 * if the distribution is truncated below 0 and renormalized
 * 
 *  a and b should both be greater than or equal to 0 but this is not checked
 */
public static double probabilityTruncZero(double a, double b, double mu, double sigma) {
    // clip at zero
    a = Math.max(a, 0.0);
    b = Math.max(b, 0.0);
    final double denom = sigma * SQRT2;
    final double scaledSDA = (a - mu) / denom; 
    final double scaledSDB = (b - mu) / denom; 
    // compute prob 
    final double probNormTimes2 = Erf.erf(scaledSDA, scaledSDB); 
    // renormalize
    final double scaledSD0 = -mu / denom; 
    final double reZTimes2 = Erf.erfc(scaledSD0);
    return probNormTimes2 / reZTimes2; 
}
 

/**
 * Merge a scan into a merged spectrum.
 * 
 * @param orderedByMz peaks from merged spectrum, sorted by ascending m/z
 * @param orderedByInt peaks from scan, sorted by descending intensity
 * @return a merged spectrum. Might be the original one if no new peaks were added.
 */
private static MergedDataPoint[] merge(MergedDataPoint[] orderedByMz, DataPoint[] orderedByInt,
    MzMergeMode mzMergeMode, IntensityMergeMode intensityMergeMode, MZTolerance expectedPPM) {
  // we assume a rather large deviation as signal peaks should be
  // contained in more than one
  // measurement
  expectedPPM =
      new MZTolerance(expectedPPM.getMzTolerance() * 4, expectedPPM.getPpmTolerance() * 4);
  final List<MergedDataPoint> append = new ArrayList<>();
  for (int k = 0; k < orderedByInt.length; ++k) {
    final DataPoint peak = orderedByInt[k];
    final double dev = expectedPPM.getMzToleranceForMass(peak.getMZ());
    final double lb = peak.getMZ() - dev, ub = peak.getMZ() + dev;
    int mz1 =
        Arrays.binarySearch(orderedByMz, peak, Comparator.comparingDouble(DataPoint::getMZ));
    if (mz1 < 0) {
      mz1 = -(mz1 + 1);
    }
    int mz0 = mz1 - 1;
    while (mz1 < orderedByMz.length && orderedByMz[mz1].getMZ() <= ub)
      ++mz1;
    --mz1;
    while (mz0 >= 0 && orderedByMz[mz0].getMZ() >= lb)
      --mz0;
    ++mz0;
    if (mz0 <= mz1) {
      // merge!
      int mostIntense = mz0;
      double bestScore = Double.NEGATIVE_INFINITY;
      for (int i = mz0; i <= mz1; ++i) {
        final double massDiff = orderedByMz[i].getMZ() - peak.getMZ();
        final double score =
            Erf.erfc(3 * massDiff) / (dev * Math.sqrt(2)) * orderedByMz[i].getIntensity();
        if (score > bestScore) {
          bestScore = score;
          mostIntense = i;
        }
      }

      orderedByMz[mostIntense] =
          orderedByMz[mostIntense].merge(peak, mzMergeMode, intensityMergeMode);

    } else {
      // append
      append.add(new MergedDataPoint(mzMergeMode, intensityMergeMode, peak));
    }
  }
  if (append.size() > 0) {
    int offset = orderedByMz.length;
    orderedByMz = Arrays.copyOf(orderedByMz, orderedByMz.length + append.size());
    for (MergedDataPoint p : append) {
      orderedByMz[offset++] = p;
    }
    ScanUtils.sortDataPointsByMz(orderedByMz);
  }
  return orderedByMz;
}
 

/**
 * Merge a scan into a merged spectrum.
 * 
 * @param orderedByMz peaks from merged spectrum, sorted by ascending m/z
 * @param orderedByInt peaks from scan, sorted by descending intensity
 * @return a merged spectrum. Might be the original one if no new peaks were added.
 */
private static MergedDataPoint[] merge(MergedDataPoint[] orderedByMz, DataPoint[] orderedByInt,
    MzMergeMode mzMergeMode, IntensityMergeMode intensityMergeMode, MZTolerance expectedPPM) {
  // we assume a rather large deviation as signal peaks should be contained in more than one
  // measurement
  expectedPPM =
      new MZTolerance(expectedPPM.getMzTolerance() * 4, expectedPPM.getPpmTolerance() * 4);
  final List<MergedDataPoint> append = new ArrayList<>();
  for (int k = 0; k < orderedByInt.length; ++k) {
    final DataPoint peak = orderedByInt[k];
    final double dev = expectedPPM.getMzToleranceForMass(peak.getMZ());
    final double lb = peak.getMZ() - dev, ub = peak.getMZ() + dev;
    int mz1 =
        Arrays.binarySearch(orderedByMz, peak, Comparator.comparingDouble(DataPoint::getMZ));
    if (mz1 < 0) {
      mz1 = -(mz1 + 1);
    }
    int mz0 = mz1 - 1;
    while (mz1 < orderedByMz.length && orderedByMz[mz1].getMZ() <= ub)
      ++mz1;
    --mz1;
    while (mz0 >= 0 && orderedByMz[mz0].getMZ() >= lb)
      --mz0;
    ++mz0;
    if (mz0 <= mz1) {
      // merge!
      int mostIntense = mz0;
      double bestScore = Double.NEGATIVE_INFINITY;
      for (int i = mz0; i <= mz1; ++i) {
        final double massDiff = orderedByMz[i].getMZ() - peak.getMZ();
        final double score =
            Erf.erfc(3 * massDiff) / (dev * Math.sqrt(2)) * orderedByMz[i].getIntensity();
        if (score > bestScore) {
          bestScore = score;
          mostIntense = i;
        }
      }

      orderedByMz[mostIntense] =
          orderedByMz[mostIntense].merge(peak, mzMergeMode, intensityMergeMode);

    } else {
      // append
      append.add(new MergedDataPoint(mzMergeMode, intensityMergeMode, peak));
    }
  }
  if (append.size() > 0) {
    int offset = orderedByMz.length;
    orderedByMz = Arrays.copyOf(orderedByMz, orderedByMz.length + append.size());
    for (MergedDataPoint p : append) {
      orderedByMz[offset++] = p;
    }
    ScanUtils.sortDataPointsByMz(orderedByMz);
  }
  return orderedByMz;
}
 

public static double cumulativeNonTrunc(double x, double mu, double sigma) {
    final double denom = sigma * SQRT2;
    return 0.5 * Erf.erfc((mu - x) / denom);
}
 

/** {@inheritDoc}
 * <p>
 * From Wikipedia: the cumulative distribution function is
 * </p>
 * <pre>
 * f(x; u, c) = erfc (&radic; (c / 2 (x - u )))
 * </pre>
 */
public double cumulativeProbability(final double x) {
    if (x < mu) {
        return Double.NaN;
    }
    return Erf.erfc(FastMath.sqrt(halfC / (x - mu)));
}
 

/** {@inheritDoc}
 * <p>
 * From Wikipedia: the cumulative distribution function is
 * </p>
 * <pre>
 * f(x; u, c) = erfc (&radic; (c / 2 (x - u )))
 * </pre>
 */
public double cumulativeProbability(final double x) {
    if (x < mu) {
        return Double.NaN;
    }
    return Erf.erfc(FastMath.sqrt(halfC / (x - mu)));
}
 

/** {@inheritDoc}
 * <p>
 * From Wikipedia: the cumulative distribution function is
 * </p>
 * <pre>
 * f(x; u, c) = erfc (&radic; (c / 2 (x - u )))
 * </pre>
 */
public double cumulativeProbability(final double x) {
    if (x < mu) {
        return Double.NaN;
    }
    return Erf.erfc(FastMath.sqrt(halfC / (x - mu)));
}
 

/** {@inheritDoc}
 * <p>
 * From Wikipedia: the cumulative distribution function is
 * </p>
 * <pre>
 * f(x; u, c) = erfc (&radic; (c / 2 (x - u )))
 * </pre>
 */
public double cumulativeProbability(final double x) {
    if (x < mu) {
        return Double.NaN;
    }
    return Erf.erfc(FastMath.sqrt(halfC / (x - mu)));
}
 

/** {@inheritDoc}
 * <p>
 * From Wikipedia: the cumulative distribution function is
 * </p>
 * <pre>
 * f(x; u, c) = erfc (&radic; (c / 2 (x - u )))
 * </pre>
 */
public double cumulativeProbability(final double x) {
    if (x < mu) {
        return Double.NaN;
    }
    return Erf.erfc(FastMath.sqrt(halfC / (x - mu)));
}