List of usage examples for org.apache.commons.math3.optim.nonlinear.vector.jacobian LevenbergMarquardtOptimizer LevenbergMarquardtOptimizer
public LevenbergMarquardtOptimizer()
From source file:net.liuxuan.temp.mathtest.java
public static void main(String[] args) { final CurveFitter fitter = new CurveFitter(new LevenbergMarquardtOptimizer()); fitter.addObservedPoint(-1.00, 2.021170021833143); fitter.addObservedPoint(-0.99, 2.221135431136975); fitter.addObservedPoint(-0.98, 2.09985277659314); fitter.addObservedPoint(-0.97, 2.0211192647627025); // ... Lots of lines omitted ... fitter.addObservedPoint(0.99, -2.4345814727089854); // The degree of the polynomial is deduced from the length of the array containing // the initial guess for the coefficients of the polynomial. final double[] init = { 12.9, -3.4, 2.1 }; // 12.9 - 3.4 x + 2.1 x^2 // Compute optimal coefficients. final double[] best = fitter.fit(new PolynomialFunction.Parametric(), init); // Construct the polynomial that best fits the data. final PolynomialFunction fitted = new PolynomialFunction(best); System.out.println(fitted.value(-0.995)); ;/* w w w . ja v a 2s . co m*/ System.out.println(fitted.value(0.995)); ; System.out.println(fitted.toString()); ; System.out.println("============================================================="); PolynomialCurveFitter pcf = PolynomialCurveFitter.create(3); WeightedObservedPoints s; List<WeightedObservedPoint> points = new ArrayList<WeightedObservedPoint>(); points.add(new WeightedObservedPoint(1, -1.00, 2.021170021833143)); points.add(new WeightedObservedPoint(1, -0.99, 2.221135431136975)); points.add(new WeightedObservedPoint(1, -0.98, 2.09985277659314)); points.add(new WeightedObservedPoint(1, -0.97, 2.0211192647627025)); points.add(new WeightedObservedPoint(1, 0.99, 2.4345814727089854)); double a[] = pcf.fit(points); for (int i = 0; i < a.length; i++) { double d = a[i]; System.out.println(d); } System.out.println(compute(a, -0.995)); System.out.println(compute(a, 0.99)); System.out.println(compute(a, 0.995)); }
From source file:com.gedaeusp.domain.NonlinearCurveFitter.java
public double[] doExponentialFitWithFixedV0(double[] v, double[] t, double v0) { CurveFitter<ExponentialParametric> fitter = new CurveFitter<ExponentialParametric>( new LevenbergMarquardtOptimizer()); addObservedPointsToFitter(v, t, fitter); return fitter.fit(new ExponentialParametric(true), new double[] { v0, 1, 1 }); }
From source file:com.gedaeusp.domain.NonlinearCurveFitter.java
public double[] doBiexponentialFitWithFixedV0(double[] v, double[] t, double[] init) { CurveFitter<BiexponentialParametric> fitter = new CurveFitter<BiexponentialParametric>( new LevenbergMarquardtOptimizer()); addObservedPointsToFitter(v, t, fitter); /* Queremos garantir que a componente rpida est associada a tau1 e a1: */ double[] best = fitter.fit(new BiexponentialParametric(true), init); if (best[Biexponential.PARAM_tau1] > best[Biexponential.PARAM_tau2]) { double tmp_tau1 = best[Biexponential.PARAM_tau2]; best[Biexponential.PARAM_tau2] = best[Biexponential.PARAM_tau1]; best[Biexponential.PARAM_tau1] = tmp_tau1; double tmp_a1 = best[Biexponential.PARAM_a2]; best[Biexponential.PARAM_a2] = best[Biexponential.PARAM_a1]; best[Biexponential.PARAM_a1] = tmp_a1; }/*from w ww . j av a2 s . co m*/ return best; }
From source file:com.gedaeusp.domain.NonlinearCurveFitter.java
public double[] doExponentialFit(double[] v, double[] t) { CurveFitter<ExponentialParametric> fitter = new CurveFitter<ExponentialParametric>( new LevenbergMarquardtOptimizer()); addObservedPointsToFitter(v, t, fitter); return fitter.fit(new ExponentialParametric(), new double[] { 0, 1, 1 }); }
From source file:com.gedaeusp.domain.NonlinearCurveFitter.java
public double[] doBiexponentialFit(double[] v, double[] t, double[] init) { CurveFitter<BiexponentialParametric> fitter = new CurveFitter<BiexponentialParametric>( new LevenbergMarquardtOptimizer()); addObservedPointsToFitter(v, t, fitter); return fitter.fit(new BiexponentialParametric(), init); }
From source file:edu.uchc.octane.CalibrationDialogAstigmatism.java
double[] parabolicFit(double[] sigma) { PolynomialFitter fitter = new PolynomialFitter(new LevenbergMarquardtOptimizer()); for (int i = 0; i < sigma.length; i++) { fitter.addObservedPoint((i - sigma.length / 2) * sliceSpacing_, sigma[i]); }//from ww w . j a v a 2 s.com double[] guess = new double[3]; guess[0] = sigma[(int) (sigma.length / 2)]; guess[2] = (sigma[0] - guess[0]) / (sliceSpacing_ * sliceSpacing_ * sigma.length * sigma.length / 4); double[] results = fitter.fit(guess); return results; }
From source file:gdsc.smlm.ij.plugins.DiffusionRateTest.java
public void run(String arg) { if (!showDialog()) return;//from www . ja v a2s . com int totalSteps = settings.seconds * settings.stepsPerSecond; final double conversionFactor = 1000000.0 / (settings.pixelPitch * settings.pixelPitch); // Diffusion rate is um^2/sec. Convert to pixels per simulation frame. final double diffusionRateInPixelsPerSecond = settings.diffusionRate * conversionFactor; final double diffusionRateInPixelsPerStep = diffusionRateInPixelsPerSecond / settings.stepsPerSecond; Utils.log(TITLE + " : D = %s um^2/sec", Utils.rounded(settings.diffusionRate, 4)); Utils.log("Mean-displacement per dimension = %s nm/sec", Utils.rounded(1e3 * ImageModel.getRandomMoveDistance(settings.diffusionRate), 4)); // Convert diffusion co-efficient into the standard deviation for the random walk final double diffusionSigma = ImageModel.getRandomMoveDistance(diffusionRateInPixelsPerStep); Utils.log("Simulation step-size = %s nm", Utils.rounded(settings.pixelPitch * diffusionSigma, 4)); // Move the molecules and get the diffusion rate IJ.showStatus("Simulating ..."); final long start = System.nanoTime(); RandomGenerator random = new Well19937c(System.currentTimeMillis() + System.identityHashCode(this)); Statistics[] stats2D = new Statistics[totalSteps]; Statistics[] stats3D = new Statistics[totalSteps]; for (int j = 0; j < totalSteps; j++) { stats2D[j] = new Statistics(); stats3D[j] = new Statistics(); } SphericalDistribution dist = new SphericalDistribution(settings.confinementRadius / settings.pixelPitch); Statistics asymptote = new Statistics(); for (int i = 0; i < settings.particles; i++) { if (i % 16 == 0) IJ.showProgress(i, settings.particles); MoleculeModel m = new MoleculeModel(i, new double[3]); if (useConfinement) { // Note: When using confinement the average displacement should asymptote // at the average distance of a point from the centre of a ball. This is 3r/4. // See: http://answers.yahoo.com/question/index?qid=20090131162630AAMTUfM // The equivalent in 2D is 2r/3. However although we are plotting 2D distance // this is a projection of the 3D position onto the plane and so the particles // will not be evenly spread (there will be clustering at centre caused by the // poles) for (int j = 0; j < totalSteps; j++) { double[] xyz = m.getCoordinates(); double[] originalXyz = Arrays.copyOf(xyz, 3); for (int n = confinementAttempts; n-- > 0;) { if (settings.useGridWalk) m.walk(diffusionSigma, random); else m.move(diffusionSigma, random); if (!dist.isWithin(m.getCoordinates())) { // Reset position for (int k = 0; k < 3; k++) xyz[k] = originalXyz[k]; } else { // The move was allowed break; } } stats2D[j].add(squared(m.getCoordinates())); stats3D[j].add(distance2(m.getCoordinates())); } asymptote.add(distance(m.getCoordinates())); } else { if (settings.useGridWalk) { for (int j = 0; j < totalSteps; j++) { m.walk(diffusionSigma, random); stats2D[j].add(squared(m.getCoordinates())); stats3D[j].add(distance2(m.getCoordinates())); } } else { for (int j = 0; j < totalSteps; j++) { m.move(diffusionSigma, random); stats2D[j].add(squared(m.getCoordinates())); stats3D[j].add(distance2(m.getCoordinates())); } } } // Debug: record all the particles so they can be analysed // System.out.printf("%f %f %f\n", m.getX(), m.getY(), m.getZ()); } final double time = (System.nanoTime() - start) / 1000000.0; IJ.showStatus("Analysing results ..."); IJ.showProgress(1); if (showDiffusionExample) { showExample(totalSteps, diffusionSigma, random); } // Plot a graph of mean squared distance double[] xValues = new double[stats2D.length]; double[] yValues2D = new double[stats2D.length]; double[] yValues3D = new double[stats3D.length]; double[] upper = new double[stats2D.length]; double[] lower = new double[stats2D.length]; final CurveFitter<Parametric> fitter2D = new CurveFitter<Parametric>(new LevenbergMarquardtOptimizer()); final CurveFitter<Parametric> fitter3D = new CurveFitter<Parametric>(new LevenbergMarquardtOptimizer()); Statistics gradient2D = new Statistics(); Statistics gradient3D = new Statistics(); final int firstN = (useConfinement) ? fitN : totalSteps; for (int j = 0; j < totalSteps; j++) { // Convert steps to seconds xValues[j] = (double) (j + 1) / settings.stepsPerSecond; // Convert values in pixels^2 to um^2 final double mean = stats2D[j].getMean() / conversionFactor; final double sd = stats2D[j].getStandardDeviation() / conversionFactor; yValues2D[j] = mean; yValues3D[j] = stats3D[j].getMean() / conversionFactor; upper[j] = mean + sd; lower[j] = mean - sd; if (j < firstN) { fitter2D.addObservedPoint(xValues[j], yValues2D[j]); gradient2D.add(yValues2D[j] / xValues[j]); fitter3D.addObservedPoint(xValues[j], yValues3D[j]); gradient3D.add(yValues3D[j] / xValues[j]); } } // TODO - Fit using the equation for 2D confined diffusion: // MSD = 4s^2 + R^2 (1 - 0.99e^(-1.84^2 Dt / R^2) // s = localisation precision // R = confinement radius // D = 2D diffusion coefficient // t = time // Do linear regression to get diffusion rate final double[] init2D = { 0, 1 / gradient2D.getMean() }; // a - b x final double[] best2D = fitter2D.fit(new PolynomialFunction.Parametric(), init2D); final PolynomialFunction fitted2D = new PolynomialFunction(best2D); final double[] init3D = { 0, 1 / gradient3D.getMean() }; // a - b x final double[] best3D = fitter3D.fit(new PolynomialFunction.Parametric(), init3D); final PolynomialFunction fitted3D = new PolynomialFunction(best3D); // Create plot String title = TITLE; Plot2 plot = new Plot2(title, "Time (seconds)", "Mean-squared Distance (um^2)", xValues, yValues2D); double[] limits = Maths.limits(upper); limits = Maths.limits(limits, lower); limits = Maths.limits(limits, yValues3D); plot.setLimits(0, totalSteps / settings.stepsPerSecond, limits[0], limits[1]); plot.setColor(Color.blue); plot.addPoints(xValues, lower, Plot2.LINE); plot.addPoints(xValues, upper, Plot2.LINE); plot.setColor(Color.magenta); plot.addPoints(xValues, yValues3D, Plot2.LINE); plot.setColor(Color.red); plot.addPoints(new double[] { xValues[0], xValues[xValues.length - 1] }, new double[] { fitted2D.value(xValues[0]), fitted2D.value(xValues[xValues.length - 1]) }, Plot2.LINE); plot.setColor(Color.green); plot.addPoints(new double[] { xValues[0], xValues[xValues.length - 1] }, new double[] { fitted3D.value(xValues[0]), fitted3D.value(xValues[xValues.length - 1]) }, Plot2.LINE); plot.setColor(Color.black); Utils.display(title, plot); // For 2D diffusion: d^2 = 4D // where: d^2 = mean-square displacement double D = best2D[1] / 4.0; String msg = "2D Diffusion rate = " + Utils.rounded(D, 4) + " um^2 / sec (" + Utils.timeToString(time) + ")"; IJ.showStatus(msg); Utils.log(msg); D = best3D[1] / 6.0; Utils.log("3D Diffusion rate = " + Utils.rounded(D, 4) + " um^2 / sec (" + Utils.timeToString(time) + ")"); if (useConfinement) Utils.log("3D asymptote distance = %s nm (expected %.2f)", Utils.rounded(asymptote.getMean() * settings.pixelPitch, 4), 3 * settings.confinementRadius / 4); }
From source file:edu.ucsc.barrel.cdf_gen.SpectrumExtract.java
private static float find511(double[] x, double[] y) { GaussianFitter fitter = new GaussianFitter(new LevenbergMarquardtOptimizer()); double[] fit_params = { 10, Constants.DOUBLE_FILL, 1 }, curve = new double[PEAK_511_WIDTH - 4]; int[] high_area = new int[PEAK_511_WIDTH]; double m, b, slope = 0, this_low = 0, last_low = 0; int apex = 0, high_cnt = 0; // guess at a linear background m = (y[PEAK_511_WIDTH - 1] - y[0]) / (x[PEAK_511_WIDTH - 1] - x[0]); b = y[0] - m * x[0];/*from ww w. j a v a 2s .c o m*/ //convert y to cnts/bin_width for (int bin_i = 0; bin_i < x.length; bin_i++) { y[bin_i] /= (RAW_EDGES[2][bin_i + PEAK_511_START + 1] - RAW_EDGES[2][bin_i + PEAK_511_START]); } //take the second derivitave to find peak for (int bin_i = 2; bin_i < x.length - 2; bin_i++) { curve[bin_i - 2] = y[bin_i + 2] - (2 * y[bin_i]) + y[bin_i - 2]; } //find low point of second derivitave using moving average this_low = (curve[0] + curve[1] + curve[2]); last_low = this_low; for (int bin_i = 2; bin_i < curve.length - 1; bin_i++) { this_low += (curve[bin_i + 1] - curve[bin_i - 2]); if (this_low < last_low) { apex = bin_i + 2; last_low = this_low; } } //do the curve fit try { fit_params[1] = x[apex]; //guess for peak location for (int bin_i = apex - 3; bin_i < apex + 3; bin_i++) { fitter.addObservedPoint(x[bin_i], y[bin_i]); } fit_params = fitter.fit(fit_params); } catch (ArrayIndexOutOfBoundsException ex) { System.out.println( "Payload ID: " + CDF_Gen.getSetting("currentPayload") + " Date: " + CDF_Gen.getSetting("date")); System.out.println("Gaussian out of bounds: " + apex); fit_params[1] = Constants.DOUBLE_FILL; } return (float) fit_params[1]; }
From source file:gdsc.smlm.ij.plugins.MeanVarianceTest.java
public void run(String arg) { if (Utils.isExtraOptions()) { ImagePlus imp = WindowManager.getCurrentImage(); if (imp.getStackSize() > 1) { GenericDialog gd = new GenericDialog(TITLE); gd.addMessage("Perform single image analysis on the current image?"); gd.addNumericField("Bias", _bias, 0); gd.showDialog();// w w w . j ava2 s . c om if (gd.wasCanceled()) return; singleImage = true; _bias = Math.abs(gd.getNextNumber()); } else { IJ.error(TITLE, "Single-image mode requires a stack"); return; } } List<ImageSample> images; String inputDirectory = ""; if (singleImage) { IJ.showStatus("Loading images..."); images = getImages(); if (images.size() == 0) { IJ.error(TITLE, "Not enough images for analysis"); return; } } else { inputDirectory = IJ.getDirectory("Select image series ..."); if (inputDirectory == null) return; SeriesOpener series = new SeriesOpener(inputDirectory, false); series.setVariableSize(true); if (series.getNumberOfImages() < 3) { IJ.error(TITLE, "Not enough images in the selected directory"); return; } if (!IJ.showMessageWithCancel(TITLE, String.format("Analyse %d images, first image:\n%s", series.getNumberOfImages(), series.getImageList()[0]))) { return; } IJ.showStatus("Loading images"); images = getImages(series); if (images.size() < 3) { IJ.error(TITLE, "Not enough images for analysis"); return; } if (images.get(0).exposure != 0) { IJ.error(TITLE, "First image in series must have exposure 0 (Bias image)"); return; } } boolean emMode = (arg != null && arg.contains("em")); if (emMode) { // Ask the user for the camera gain ... GenericDialog gd = new GenericDialog(TITLE); gd.addMessage("Estimating the EM-gain requires the camera gain without EM readout enabled"); gd.addNumericField("Camera_gain (ADU/e-)", cameraGain, 4); gd.showDialog(); if (gd.wasCanceled()) return; cameraGain = gd.getNextNumber(); } IJ.showStatus("Computing mean & variance"); final double nImages = images.size(); for (int i = 0; i < images.size(); i++) { IJ.showStatus(String.format("Computing mean & variance %d/%d", i + 1, images.size())); images.get(i).compute(singleImage, i / nImages, (i + 1) / nImages); } IJ.showProgress(1); IJ.showStatus("Computing results"); TextWindow results = createResultsWindow(); // Allow user to input multiple bias images int start = 0; Statistics biasStats = new Statistics(); Statistics noiseStats = new Statistics(); final double bias; if (singleImage) { bias = _bias; } else { while (start < images.size()) { ImageSample sample = images.get(start); if (sample.exposure == 0) { biasStats.add(sample.means); for (PairSample pair : sample.samples) { noiseStats.add(pair.variance); } start++; } else break; } bias = biasStats.getMean(); } // Get the mean-variance data int total = 0; for (int i = start; i < images.size(); i++) total += images.get(i).samples.size(); double[] mean = new double[total]; double[] variance = new double[mean.length]; Statistics gainStats = (singleImage) ? new StoredDataStatistics(total) : new Statistics(); final CurveFitter<Parametric> fitter = new CurveFitter<Parametric>(new LevenbergMarquardtOptimizer()); for (int i = (singleImage) ? 0 : start, j = 0; i < images.size(); i++) { ImageSample sample = images.get(i); for (PairSample pair : sample.samples) { if (j % 16 == 0) IJ.showProgress(j, total); mean[j] = pair.getMean(); variance[j] = pair.variance; // Gain is in ADU / e double gain = variance[j] / (mean[j] - bias); gainStats.add(gain); fitter.addObservedPoint(mean[j], variance[j]); if (emMode) { gain /= (2 * cameraGain); } StringBuilder sb = new StringBuilder(); sb.append(sample.title).append("\t"); sb.append(sample.exposure).append("\t"); sb.append(pair.slice1).append("\t"); sb.append(pair.slice2).append("\t"); sb.append(IJ.d2s(pair.mean1, 2)).append("\t"); sb.append(IJ.d2s(pair.mean2, 2)).append("\t"); sb.append(IJ.d2s(mean[j], 2)).append("\t"); sb.append(IJ.d2s(variance[j], 2)).append("\t"); sb.append(Utils.rounded(gain, 4)); results.append(sb.toString()); j++; } } IJ.showProgress(1); if (singleImage) { StoredDataStatistics stats = (StoredDataStatistics) gainStats; Utils.log(TITLE); if (emMode) { double[] values = stats.getValues(); MathArrays.scaleInPlace(0.5, values); stats = new StoredDataStatistics(values); } // Plot the gain over time String title = TITLE + " Gain vs Frame"; Plot2 plot = new Plot2(title, "Slice", "Gain", Utils.newArray(gainStats.getN(), 1, 1.0), stats.getValues()); PlotWindow pw = Utils.display(title, plot); // Show a histogram String label = String.format("Mean = %s, Median = %s", Utils.rounded(stats.getMean()), Utils.rounded(stats.getMedian())); int id = Utils.showHistogram(TITLE, stats, "Gain", 0, 1, 100, true, label); if (Utils.isNewWindow()) { Point point = pw.getLocation(); point.x = pw.getLocation().x; point.y += pw.getHeight(); WindowManager.getImage(id).getWindow().setLocation(point); } Utils.log("Single-image mode: %s camera", (emMode) ? "EM-CCD" : "Standard"); final double gain = stats.getMedian(); if (emMode) { final double totalGain = gain; final double emGain = totalGain / cameraGain; Utils.log(" Gain = 1 / %s (ADU/e-)", Utils.rounded(cameraGain, 4)); Utils.log(" EM-Gain = %s", Utils.rounded(emGain, 4)); Utils.log(" Total Gain = %s (ADU/e-)", Utils.rounded(totalGain, 4)); } else { cameraGain = gain; Utils.log(" Gain = 1 / %s (ADU/e-)", Utils.rounded(cameraGain, 4)); } } else { IJ.showStatus("Computing fit"); // Sort int[] indices = rank(mean); mean = reorder(mean, indices); variance = reorder(variance, indices); // Compute optimal coefficients. final double[] init = { 0, 1 / gainStats.getMean() }; // a - b x final double[] best = fitter.fit(new PolynomialFunction.Parametric(), init); // Construct the polynomial that best fits the data. final PolynomialFunction fitted = new PolynomialFunction(best); // Plot mean verses variance. Gradient is gain in ADU/e. String title = TITLE + " results"; Plot2 plot = new Plot2(title, "Mean", "Variance"); double[] xlimits = Maths.limits(mean); double[] ylimits = Maths.limits(variance); double xrange = (xlimits[1] - xlimits[0]) * 0.05; if (xrange == 0) xrange = 0.05; double yrange = (ylimits[1] - ylimits[0]) * 0.05; if (yrange == 0) yrange = 0.05; plot.setLimits(xlimits[0] - xrange, xlimits[1] + xrange, ylimits[0] - yrange, ylimits[1] + yrange); plot.setColor(Color.blue); plot.addPoints(mean, variance, Plot2.CROSS); plot.setColor(Color.red); plot.addPoints(new double[] { mean[0], mean[mean.length - 1] }, new double[] { fitted.value(mean[0]), fitted.value(mean[mean.length - 1]) }, Plot2.LINE); Utils.display(title, plot); final double avBiasNoise = Math.sqrt(noiseStats.getMean()); Utils.log(TITLE); Utils.log(" Directory = %s", inputDirectory); Utils.log(" Bias = %s +/- %s (ADU)", Utils.rounded(bias, 4), Utils.rounded(avBiasNoise, 4)); Utils.log(" Variance = %s + %s * mean", Utils.rounded(best[0], 4), Utils.rounded(best[1], 4)); if (emMode) { final double emGain = best[1] / (2 * cameraGain); // Noise is standard deviation of the bias image divided by the total gain (in ADU/e-) final double totalGain = emGain * cameraGain; Utils.log(" Read Noise = %s (e-) [%s (ADU)]", Utils.rounded(avBiasNoise / totalGain, 4), Utils.rounded(avBiasNoise, 4)); Utils.log(" Gain = 1 / %s (ADU/e-)", Utils.rounded(1 / cameraGain, 4)); Utils.log(" EM-Gain = %s", Utils.rounded(emGain, 4)); Utils.log(" Total Gain = %s (ADU/e-)", Utils.rounded(totalGain, 4)); } else { // Noise is standard deviation of the bias image divided by the gain (in ADU/e-) cameraGain = best[1]; final double readNoise = avBiasNoise / cameraGain; Utils.log(" Read Noise = %s (e-) [%s (ADU)]", Utils.rounded(readNoise, 4), Utils.rounded(readNoise * cameraGain, 4)); Utils.log(" Gain = 1 / %s (ADU/e-)", Utils.rounded(1 / cameraGain, 4)); } } IJ.showStatus(""); }
From source file:gdsc.smlm.fitting.BinomialFitter.java
/** * Fit the binomial distribution (n,p) to the cumulative histogram. Performs fitting assuming a fixed n value and * attempts to optimise p./*from w w w . jav a 2s. c o m*/ * * @param histogram * The input histogram * @param mean * The histogram mean (used to estimate p). Calculated if NaN. * @param n * The n to evaluate * @param zeroTruncated * True if the model should ignore n=0 (zero-truncated binomial) * @return The best fit (n, p) * @throws IllegalArgumentException * If any of the input data values are negative * @throws IllegalArgumentException * If any fitting a zero truncated binomial and there are no values above zero */ public PointValuePair fitBinomial(double[] histogram, double mean, int n, boolean zeroTruncated) { if (Double.isNaN(mean)) mean = getMean(histogram); if (zeroTruncated && histogram[0] > 0) { log("Fitting zero-truncated histogram but there are zero values - Renormalising to ignore zero"); double cumul = 0; for (int i = 1; i < histogram.length; i++) cumul += histogram[i]; if (cumul == 0) throw new IllegalArgumentException( "Fitting zero-truncated histogram but there are no non-zero values"); histogram[0] = 0; for (int i = 1; i < histogram.length; i++) histogram[i] /= cumul; } int nFittedPoints = Math.min(histogram.length, n + 1) - ((zeroTruncated) ? 1 : 0); if (nFittedPoints < 1) { log("No points to fit (%d): Histogram.length = %d, n = %d, zero-truncated = %b", nFittedPoints, histogram.length, n, zeroTruncated); return null; } // The model is only fitting the probability p // For a binomial n*p = mean => p = mean/n double[] initialSolution = new double[] { FastMath.min(mean / n, 1) }; // Create the function BinomialModelFunction function = new BinomialModelFunction(histogram, n, zeroTruncated); double[] lB = new double[1]; double[] uB = new double[] { 1 }; SimpleBounds bounds = new SimpleBounds(lB, uB); // Fit // CMAESOptimizer or BOBYQAOptimizer support bounds // CMAESOptimiser based on Matlab code: // https://www.lri.fr/~hansen/cmaes.m // Take the defaults from the Matlab documentation int maxIterations = 2000; double stopFitness = 0; //Double.NEGATIVE_INFINITY; boolean isActiveCMA = true; int diagonalOnly = 0; int checkFeasableCount = 1; RandomGenerator random = new Well19937c(); boolean generateStatistics = false; ConvergenceChecker<PointValuePair> checker = new SimpleValueChecker(1e-6, 1e-10); // The sigma determines the search range for the variables. It should be 1/3 of the initial search region. OptimizationData sigma = new CMAESOptimizer.Sigma(new double[] { (uB[0] - lB[0]) / 3 }); OptimizationData popSize = new CMAESOptimizer.PopulationSize((int) (4 + Math.floor(3 * Math.log(2)))); try { PointValuePair solution = null; boolean noRefit = maximumLikelihood; if (n == 1 && zeroTruncated) { // No need to fit solution = new PointValuePair(new double[] { 1 }, 0); noRefit = true; } else { GoalType goalType = (maximumLikelihood) ? GoalType.MAXIMIZE : GoalType.MINIMIZE; // Iteratively fit CMAESOptimizer opt = new CMAESOptimizer(maxIterations, stopFitness, isActiveCMA, diagonalOnly, checkFeasableCount, random, generateStatistics, checker); for (int iteration = 0; iteration <= fitRestarts; iteration++) { try { // Start from the initial solution PointValuePair result = opt.optimize(new InitialGuess(initialSolution), new ObjectiveFunction(function), goalType, bounds, sigma, popSize, new MaxIter(maxIterations), new MaxEval(maxIterations * 2)); //System.out.printf("CMAES Iter %d initial = %g (%d)\n", iteration, result.getValue(), // opt.getEvaluations()); if (solution == null || result.getValue() < solution.getValue()) { solution = result; } } catch (TooManyEvaluationsException e) { } catch (TooManyIterationsException e) { } if (solution == null) continue; try { // Also restart from the current optimum PointValuePair result = opt.optimize(new InitialGuess(solution.getPointRef()), new ObjectiveFunction(function), goalType, bounds, sigma, popSize, new MaxIter(maxIterations), new MaxEval(maxIterations * 2)); //System.out.printf("CMAES Iter %d restart = %g (%d)\n", iteration, result.getValue(), // opt.getEvaluations()); if (result.getValue() < solution.getValue()) { solution = result; } } catch (TooManyEvaluationsException e) { } catch (TooManyIterationsException e) { } } if (solution == null) return null; } if (noRefit) { // Although we fit the log-likelihood, return the sum-of-squares to allow // comparison across different n double p = solution.getPointRef()[0]; double ss = 0; double[] obs = function.p; double[] exp = function.getP(p); for (int i = 0; i < obs.length; i++) ss += (obs[i] - exp[i]) * (obs[i] - exp[i]); return new PointValuePair(solution.getPointRef(), ss); } // We can do a LVM refit if the number of fitted points is more than 1 else if (nFittedPoints > 1) { // Improve SS fit with a gradient based LVM optimizer LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer(); try { final BinomialModelFunctionGradient gradientFunction = new BinomialModelFunctionGradient( histogram, n, zeroTruncated); PointVectorValuePair lvmSolution = optimizer.optimize(new MaxIter(3000), new MaxEval(Integer.MAX_VALUE), new ModelFunctionJacobian(new MultivariateMatrixFunction() { public double[][] value(double[] point) throws IllegalArgumentException { return gradientFunction.jacobian(point); } }), new ModelFunction(gradientFunction), new Target(gradientFunction.p), new Weight(gradientFunction.getWeights()), new InitialGuess(solution.getPointRef())); double ss = 0; double[] obs = gradientFunction.p; double[] exp = lvmSolution.getValue(); for (int i = 0; i < obs.length; i++) ss += (obs[i] - exp[i]) * (obs[i] - exp[i]); // Check the pValue is valid since the LVM is not bounded. double p = lvmSolution.getPointRef()[0]; if (ss < solution.getValue() && p <= 1 && p >= 0) { //log("Re-fitting improved the SS from %s to %s (-%s%%)", // Utils.rounded(solution.getValue(), 4), Utils.rounded(ss, 4), // Utils.rounded(100 * (solution.getValue() - ss) / solution.getValue(), 4)); return new PointValuePair(lvmSolution.getPoint(), ss); } } catch (TooManyIterationsException e) { log("Failed to re-fit: Too many iterations (%d)", optimizer.getIterations()); } catch (ConvergenceException e) { log("Failed to re-fit: %s", e.getMessage()); } catch (Exception e) { // Ignore this ... } } return solution; } catch (Exception e) { log("Failed to fit Binomial distribution with N=%d : %s", n, e.getMessage()); } return null; }