Example usage for org.apache.commons.math3.exception.util LocalizedFormats BANDWIDTH

List of usage examples for org.apache.commons.math3.exception.util LocalizedFormats BANDWIDTH

Introduction

In this page you can find the example usage for org.apache.commons.math3.exception.util LocalizedFormats BANDWIDTH.

Prototype

LocalizedFormats BANDWIDTH

To view the source code for org.apache.commons.math3.exception.util LocalizedFormats BANDWIDTH.

Click Source Link

Usage

From source file:cz.cuni.lf1.lge.ThunderSTORM.results.ModifiedLoess.java

/**
 * Construct a new {@link LoessInterpolator}
 * with given bandwidth, number of robustness iterations and accuracy.
 *
 * @param bandwidth  when computing the loess fit at
 * a particular point, this fraction of source points closest
 * to the current point is taken into account for computing
 * a least-squares regression.</br>
 * A sensible value is usually 0.25 to 0.5, the default value is
 * {@link #DEFAULT_BANDWIDTH}./*from w w  w  . j a  v a2s .  c  o m*/
 * @param robustnessIters This many robustness iterations are done.</br>
 * A sensible value is usually 0 (just the initial fit without any
 * robustness iterations) to 4, the default value is
 * {@link #DEFAULT_ROBUSTNESS_ITERS}.
 * @param accuracy If the median residual at a certain robustness iteration
 * is less than this amount, no more iterations are done.
 * @throws OutOfRangeException if bandwidth does not lie in the interval [0,1].
 * @throws NotPositiveException if {@code robustnessIters} is negative.
 * @see #LoessInterpolator(double, int)
 * @since 2.1
 */
public ModifiedLoess(double bandwidth, int robustnessIters, double accuracy)
        throws OutOfRangeException, NotPositiveException {
    if (bandwidth < 0 || bandwidth > 1) {
        throw new OutOfRangeException(LocalizedFormats.BANDWIDTH, bandwidth, 0, 1);
    }
    this.bandwidth = bandwidth;
    if (robustnessIters < 0) {
        throw new NotPositiveException(LocalizedFormats.ROBUSTNESS_ITERATIONS, robustnessIters);
    }
    this.robustnessIters = robustnessIters;
    this.accuracy = accuracy;
}

From source file:cz.cuni.lf1.lge.ThunderSTORM.results.ModifiedLoess.java

/**
 * Compute a weighted loess fit on the data at the original abscissae.
 *
 * @param xval Arguments for the interpolation points.
 * @param yval Values for the interpolation points.
 * @param weights point weights: coefficients by which the robustness weight
 * of a point is multiplied./*from  www .  j a va2 s . c o  m*/
 * @return the values of the loess fit at corresponding original abscissae.
 * @throws NonMonotonicSequenceException if {@code xval} not sorted in
 * strictly increasing order.
 * @throws DimensionMismatchException if {@code xval} and {@code yval} have
 * different sizes.
 * @throws NoDataException if {@code xval} or {@code yval} has zero size.
 * @throws NotFiniteNumberException if any of the arguments and values are
 not finite real numbers.
 * @throws NumberIsTooSmallException if the bandwidth is too small to
 * accomodate the size of the input data (i.e. the bandwidth must be
 * larger than 2/n).
 * @since 2.1
 */
public final double[] smooth(final double[] xval, final double[] yval, final double[] weights)
        throws NonMonotonicSequenceException, DimensionMismatchException, NoDataException,
        NotFiniteNumberException, NumberIsTooSmallException {
    if (xval.length != yval.length) {
        throw new DimensionMismatchException(xval.length, yval.length);
    }

    final int n = xval.length;

    if (n == 0) {
        throw new NoDataException();
    }

    checkAllFiniteReal(xval);
    checkAllFiniteReal(yval);
    checkAllFiniteReal(weights);

    MathArrays.checkOrder(xval, MathArrays.OrderDirection.INCREASING, false);

    if (n == 1) {
        return new double[] { yval[0] };
    }

    if (n == 2) {
        return new double[] { yval[0], yval[1] };
    }

    int bandwidthInPoints = (int) (bandwidth * n);

    if (bandwidthInPoints < 2) {
        throw new NumberIsTooSmallException(LocalizedFormats.BANDWIDTH, bandwidthInPoints, 2, true);
    }

    final double[] res = new double[n];

    final double[] residuals = new double[n];
    final double[] sortedResiduals = new double[n];

    final double[] robustnessWeights = new double[n];

    // Do an initial fit and 'robustnessIters' robustness iterations.
    // This is equivalent to doing 'robustnessIters+1' robustness iterations
    // starting with all robustness weights set to 1.
    Arrays.fill(robustnessWeights, 1);

    for (int iter = 0; iter <= robustnessIters; ++iter) {
        final int[] bandwidthInterval = { 0, bandwidthInPoints - 1 };
        // At each x, compute a local weighted linear regression
        for (int i = 0; i < n; ++i) {
            final double x = xval[i];

            // Find out the interval of source points on which
            // a regression is to be made.
            if (i > 0) {
                updateBandwidthInterval(xval, weights, i, bandwidthInterval);
            }

            final int ileft = bandwidthInterval[0];
            final int iright = bandwidthInterval[1];

            // Compute the point of the bandwidth interval that is
            // farthest from x
            final int edge;
            if (xval[i] - xval[ileft] > xval[iright] - xval[i]) {
                edge = ileft;
            } else {
                edge = iright;
            }

            // Compute a least-squares linear fit weighted by
            // the product of robustness weights and the tricube
            // weight function.
            // See http://en.wikipedia.org/wiki/Linear_regression
            // (section "Univariate linear case")
            // and http://en.wikipedia.org/wiki/Weighted_least_squares
            // (section "Weighted least squares")
            double sumWeights = 0;
            double sumX = 0;
            double sumXSquared = 0;
            double sumY = 0;
            double sumXY = 0;
            double denom = FastMath.abs(1.0 / (xval[edge] - x));
            for (int k = ileft; k <= iright; ++k) {
                final double xk = xval[k];
                final double yk = yval[k];
                final double dist = (k < i) ? x - xk : xk - x;
                final double w = tricube(dist * denom) * robustnessWeights[k] * weights[k];
                final double xkw = xk * w;
                sumWeights += w;
                sumX += xkw;
                sumXSquared += xk * xkw;
                sumY += yk * w;
                sumXY += yk * xkw;
            }

            final double meanX = sumX / sumWeights;
            final double meanY = sumY / sumWeights;
            final double meanXY = sumXY / sumWeights;
            final double meanXSquared = sumXSquared / sumWeights;

            final double beta;
            if (FastMath.sqrt(FastMath.abs(meanXSquared - meanX * meanX)) < accuracy) {
                beta = 0;
            } else {
                beta = (meanXY - meanX * meanY) / (meanXSquared - meanX * meanX);
            }

            final double alpha = meanY - beta * meanX;

            res[i] = beta * x + alpha;
            residuals[i] = FastMath.abs(yval[i] - res[i]);
        }

        // No need to recompute the robustness weights at the last
        // iteration, they won't be needed anymore
        if (iter == robustnessIters) {
            break;
        }

        // Recompute the robustness weights.

        // Find the median residual.
        // An arraycopy and a sort are completely tractable here,
        // because the preceding loop is a lot more expensive
        System.arraycopy(residuals, 0, sortedResiduals, 0, n);
        Arrays.sort(sortedResiduals);
        final double medianResidual = sortedResiduals[n / 2];

        if (FastMath.abs(medianResidual) < accuracy) {
            break;
        }

        for (int i = 0; i < n; ++i) {
            final double arg = residuals[i] / (6 * medianResidual);
            if (arg >= 1) {
                robustnessWeights[i] = 0;
            } else {
                final double w = 1 - arg * arg;
                robustnessWeights[i] = w * w;
            }
        }
    }

    return res;
}