List of usage examples for org.apache.commons.math3.util MathArrays checkOrder
public static void checkOrder(double[] val, OrderDirection dir, boolean strict) throws NonMonotonicSequenceException
From source file:gedi.util.math.function.KnotFunction.java
/** * Builds a knot function from a list of arguments and the corresponding * values. Specifically, returns the function h(x) defined by <pre><code> * h(x) = y[0] for all x < x[1]/* ww w.j a v a 2s . co m*/ * y[1] for x[1] <= x < x[2] * ... * y[y.length - 1] for x >= x[x.length - 1] * </code></pre> * The value of {@code x[0]} is ignored, but it must be strictly less than * {@code x[1]}. * * x and y are not copied! * if x and y are not strictly monotonically increasing, the mean of the y values is computed for equal x values * It has to be increasing! * * @param x Domain values where the function changes value. * @param y Values of the function. * @throws NullArgumentException if {@code x} or {@code y} are {@code null}. * @throws NoDataException if {@code x} or {@code y} are zero-length. * @throws DimensionMismatchException if {@code x} and {@code y} do not * have the same length. */ public KnotFunction(double[] x, double[] y) throws NullArgumentException, NoDataException, DimensionMismatchException { if (x == null || y == null) { throw new NullArgumentException(); } if (x.length == 0 || y.length == 0) { throw new NoDataException(); } if (y.length != x.length) { throw new DimensionMismatchException(y.length, x.length); } MathArrays.checkOrder(x, OrderDirection.INCREASING, false); int index = 0; double sumy = 0; int ny = 0; for (int i = 0; i < x.length; i++) { if (x[i] != x[index]) { index++; ny = 1; sumy = y[i]; } else { ny++; sumy += y[i]; } x[index] = x[i]; y[index] = sumy / ny; } if (index + 1 != x.length) { x = Arrays.copyOf(x, index + 1); y = Arrays.copyOf(y, index + 1); } this.x = x; this.y = y; }
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./*w w w . j a v a 2 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; }