Examples with WeightedObservedPoint org.apache.commons.math3.fitting.WeightedObservedPoint used on opensource projects

Example 6 with WeightedObservedPoint

use of org.apache.commons.math3.fitting.WeightedObservedPoint in project ffx by mjschnie.

the class BlockAverager method residual.

/**
 * Computes a residual to the given points for the provided fit type and
 * coefficients.
 */
private double residual(List<WeightedObservedPoint> obs, FITTER fitter, double[] coeffs) {
    int n = obs.size();
    double sumydiffsq = 0.0;
    for (int i = 0; i < n; i++) {
        final double x = obs.get(i).getX() * obs.get(i).getWeight();
        final double y = obs.get(i).getY() * obs.get(i).getWeight();
        double value;
        switch(fitter) {
            case LINEAR:
                value = coeffs[0] + coeffs[1] * x;
                break;
            case LOG:
                value = coeffs[0] + coeffs[1] * log(x);
                break;
            case POWER:
                value = coeffs[0] * pow(x, coeffs[1]);
                break;
            default:
                throw new UnsupportedOperationException();
        }
        sumydiffsq += (y - value) * (y - value);
    }
    double residual = sqrt(sumydiffsq);
    return residual;
}

Example 7 with WeightedObservedPoint

use of org.apache.commons.math3.fitting.WeightedObservedPoint in project ffx by mjschnie.

the class BlockAverager method computeBinUncertainties.

/**
 * Compute the statistical uncertainty of G in each histogram bin and
 * overall. Loop over increasing values of block size. For each, calculate
 * the block means and their standard deviation. Then limit(blockStdErr,
 * blockSize to entireTraj) == trajStdErr.
 *
 * @return aggregate standard error of the total free energy change
 */
public double[] computeBinUncertainties() {
    double[][] sems = new double[numBins][maxBlockSize + 1];
    BinDev[][] binStDevs = new BinDev[numBins][maxBlockSize + 1];
    List<WeightedObservedPoint>[] obsDev = new ArrayList[numBins];
    List<WeightedObservedPoint>[] obsErr = new ArrayList[numBins];
    for (int binIndex = 0; binIndex < numBins; binIndex++) {
        logger.info(format(" Computing stdError for bin %d...", binIndex));
        obsDev[binIndex] = new ArrayList<>();
        obsErr[binIndex] = new ArrayList<>();
        for (int blockSize = 1; blockSize <= maxBlockSize; blockSize += blockSizeStep) {
            int numBlocks = (int) Math.floor(numObs / blockSize);
            binStDevs[binIndex][blockSize] = new BinDev(binIndex, blockSize);
            sems[binIndex][blockSize] = binStDevs[binIndex][blockSize].stdev / Math.sqrt(numBlocks - 1);
            obsDev[binIndex].add(new WeightedObservedPoint(1.0, blockSize, binStDevs[binIndex][blockSize].stdev));
            obsErr[binIndex].add(new WeightedObservedPoint(1.0, blockSize, sems[binIndex][blockSize]));
            if (TEST) {
                logger.info(format("  bin,blockSize,stdev,sem: %d %d %.6g %.6g", binIndex, blockSize, binStDevs[binIndex][blockSize].stdev, sems[binIndex][blockSize]));
            }
        }
    }
    // Fit a function to (blockSize v. stdError) and extrapolate to blockSize == entire trajectory.
    // This is our correlation-corrected estimate of the std error for this lambda bin.
    stdError = new double[numBins];
    for (int binIndex = 0; binIndex < numBins; binIndex++) {
        logger.info(format("\n Bin %d : fitting & extrapolating blockSize v. stdError", binIndex));
        if (fitter == FITTER.POLYNOMIAL) {
            // Fit a polynomial (shitty).
            double[] coeffsPoly = PolynomialCurveFitter.create(polyDegree).fit(obsErr[binIndex]);
            PolynomialFunction poly = new PolynomialFunction(coeffsPoly);
            logger.info(format("    Poly %d:   %12.10g     %s", polyDegree, poly.value(numObs), Arrays.toString(poly.getCoefficients())));
        } else if (fitter == FITTER.POWER) {
            // Fit a power function (better).
            double[] coeffsPow = powerFit(obsErr[binIndex]);
            double powerExtrapolated = coeffsPow[0] * pow(numObs, coeffsPow[1]);
            logger.info(format("    Power:     %12.10g     %s", powerExtrapolated, Arrays.toString(coeffsPow)));
        }
        // Fit a log function (best).
        double[] logCoeffs = logFit(obsErr[binIndex]);
        double logExtrap = logCoeffs[0] + logCoeffs[1] * log(numObs);
        logger.info(format("    Log sem:   %12.10g     Residual: %12.10g     Coeffs: %6.4f, %6.4f \n", logExtrap, logCoeffs[0], logCoeffs[1]));
        // Also try fitting a linear function for the case of uncorrelated or extremely well-converged data.
        double[] linearCoef = linearFit(obsErr[binIndex]);
        double linearExtrap = linearCoef[0] + linearCoef[1] * numObs;
        logger.info(format("    Lin. sem:  %12.10g     Residual: %12.10g     Coeffs: %6.4f, %6.4f \n", linearExtrap, linearCoef[0], linearCoef[1]));
        stdError[binIndex] = logExtrap;
    }
    return stdError;
}

Example 8 with WeightedObservedPoint

use of org.apache.commons.math3.fitting.WeightedObservedPoint in project ffx by mjschnie.

the class BlockAverager method linearFit.

/**
 * Returns [A,B,R^2] such that f(x)= a + b*x is minimized for the input
 * points. As described at
 * http://mathworld.wolfram.com/LeastSquaresFitting.html
 */
private double[] linearFit(List<WeightedObservedPoint> obs) {
    int n = obs.size();
    double sumx = 0.0;
    double sumy = 0.0;
    double sumxsq = 0.0;
    double sumysq = 0.0;
    double sumxy = 0.0;
    for (int i = 0; i < n; i++) {
        final double x = obs.get(i).getX() * obs.get(i).getWeight();
        final double y = obs.get(i).getY() * obs.get(i).getWeight();
        sumx += x;
        sumy += y;
        sumxsq += x * x;
        sumysq += y * y;
        sumxy += x * y;
    }
    final double xbar = sumx / n;
    final double ybar = sumy / n;
    final double ssxx = sumxsq - n * xbar * xbar;
    final double ssyy = sumysq - n * ybar * ybar;
    final double ssxy = sumxy - n * xbar * ybar;
    final double rsq = (ssxy * ssxy) / (ssxx * ssyy);
    final double a = (ybar * sumxsq - xbar * sumxy) / (sumxsq - n * xbar * xbar);
    final double b = (sumxy - n * xbar * ybar) / (sumxsq - n * xbar * xbar);
    double[] ret = { a, b, rsq };
    return ret;
}

Example 9 with WeightedObservedPoint

use of org.apache.commons.math3.fitting.WeightedObservedPoint in project GDSC-SMLM by aherbert.

the class Fire method fitGaussian.

/**
 * Fit gaussian.
 *
 * @param x the x
 * @param y the y
 * @return new double[] { norm, mean, sigma }
 */
private static double[] fitGaussian(float[] x, float[] y) {
    final WeightedObservedPoints obs = new WeightedObservedPoints();
    for (int i = 0; i < x.length; i++) {
        obs.add(x[i], y[i]);
    }
    final Collection<WeightedObservedPoint> observations = obs.toList();
    final GaussianCurveFitter fitter = GaussianCurveFitter.create().withMaxIterations(2000);
    final GaussianCurveFitter.ParameterGuesser guess = new GaussianCurveFitter.ParameterGuesser(observations);
    double[] initialGuess = null;
    try {
        initialGuess = guess.guess();
        return fitter.withStartPoint(initialGuess).fit(observations);
    } catch (final Exception ex) {
    // We are expecting TooManyEvaluationsException.
    // Just in case there is another exception type, or the initial estimate failed
    // we catch all exceptions.
    }
    return initialGuess;
}