List of usage examples for org.apache.commons.math3.analysis UnivariateFunction UnivariateFunction
UnivariateFunction
From source file:fsm.series.Series_SS.java
public UnivariateFunction getFunction(int m) { UnivariateFunction f = new UnivariateFunction() { @Override/* w ww . jav a 2s . c o m*/ public double value(double x) { return getFunctionValue(x, m); } }; return f; }
From source file:eu.crisis_economics.abm.algorithms.statistics.FloorInterpolator.java
@Override public UnivariateFunction interpolate(final double[] x, final double[] y) throws MathIllegalArgumentException, DimensionMismatchException { Preconditions.checkNotNull(x, y);/*from w w w . ja va 2s .c om*/ if (x.length != y.length) throw new DimensionMismatchException(x.length, y.length); final int seriesLength = x.length; if (seriesLength == 0) throw new IllegalArgumentException( "FloorInterpolator.interpolate: input data is empty. Interpolation cannot continue."); final TreeMap<Double, Double> tree = new TreeMap<Double, Double>(); for (int i = 0; i < seriesLength; ++i) tree.put(x[i], y[i]); final UnivariateFunction result = new UnivariateFunction() { @Override public double value(final double t) { if (t < x[0]) return y[0]; else if (t > x[seriesLength - 1]) return y[seriesLength - 1]; else return tree.floorEntry(t).getValue(); } }; return result; }
From source file:de.unijena.bioinf.FragmentationTreeConstruction.computation.recalibration.ExpRek.java
@Override public UnivariateFunction recalibrate(MutableSpectrum<Peak> spectrum, Spectrum<Peak> referenceSpectrum) { final SimpleMutableSpectrum ref = new SimpleMutableSpectrum(referenceSpectrum); final SimpleMutableSpectrum mes = new SimpleMutableSpectrum(spectrum); final TDoubleList list = new TDoubleArrayList(spectrum.size()); // 1. allow only peaks with >2.5% intensity for recalibration for (int i = 0; i < mes.size(); ++i) { if (mes.getIntensityAt(i) <= 0.025) { ref.removePeakAt(i);/*from w w w . java 2s .co m*/ mes.removePeakAt(i); --i; } else if (mes.getIntensityAt(i) > 0.04) { list.add(mes.getMzAt(i)); } } final Deviation dev = new Deviation(4, 1e-3); final double[][] values = MzRecalibration.maxIntervalStabbing(mes, ref, new UnivariateFunction() { @Override public double value(double x) { return dev.absoluteFor(x); } }); // 2. there have to be at least 6 peaks with > 4% intensity int found = 0; for (double x : values[0]) { if (list.contains(x)) ++found; } if (found < 6) return new Identity(); final UnivariateFunction recalibration = MzRecalibration.getLinearRecalibration(values[0], values[1]); MzRecalibration.recalibrate(spectrum, recalibration); return recalibration; }
From source file:beast.math.distribution.GammaDistributionTest.java
@Test public void testPdf() { final int numberOfTests = 100; double totErr = 0; double ptotErr = 0; int np = 0;//from w w w.j a v a2s . c o m double qtotErr = 0; Random random = new Random(37); for (int i = 0; i < numberOfTests; i++) { final double mean = .01 + (3 - 0.01) * random.nextDouble(); final double var = .01 + (3 - 0.01) * random.nextDouble(); final double scale = var / mean; final double shape = mean / scale; final GammaDistribution gamma = new GammaDistribution(shape, scale); final double value = gamma.nextGamma(); final double mypdf = mypdf(value, shape, scale); final double pdf = gamma.pdf(value); if (Double.isInfinite(mypdf) && Double.isInfinite(pdf)) { continue; } assertFalse(Double.isNaN(mypdf)); assertFalse(Double.isNaN(pdf)); totErr += mypdf != 0 ? Math.abs((pdf - mypdf) / mypdf) : pdf; assertFalse("nan", Double.isNaN(totErr)); //assertEquals("" + shape + "," + scale + "," + value, mypdf,gamma.pdf(value),1e-10); final double cdf = gamma.cdf(value); UnivariateFunction f = new UnivariateFunction() { public double value(double v) { return mypdf(v, shape, scale); } }; final UnivariateIntegrator integrator = new RombergIntegrator(MachineAccuracy.SQRT_EPSILON, 1e-14, 1, 16); double x; try { x = integrator.integrate(16, f, 0.0, value); ptotErr += cdf != 0.0 ? Math.abs(x - cdf) / cdf : x; np += 1; //assertTrue("" + shape + "," + scale + "," + value + " " + Math.abs(x-cdf)/x + "> 1e-6", Math.abs(1-cdf/x) < 1e-6); //System.out.println(shape + "," + scale + " " + value); } catch (MaxCountExceededException e) { // can't integrate , skip test // System.out.println(shape + "," + scale + " skipped"); } final double q = gamma.quantile(cdf); qtotErr += q != 0 ? Math.abs(q - value) / q : value; // assertEquals("" + shape + "," + scale + "," + value + " " + Math.abs(q-value)/value, q, value, 1e-6); } //System.out.println( !Double.isNaN(totErr) ); // System.out.println(totErr); // bad test, but I can't find a good threshold that works for all individual cases assertTrue("failed " + totErr / numberOfTests, totErr / numberOfTests < 1e-7); assertTrue("failed " + ptotErr / np, np > 0 ? (ptotErr / np < 1e-5) : true); assertTrue("failed " + qtotErr / numberOfTests, qtotErr / numberOfTests < 1e-7); }
From source file:eu.crisis_economics.abm.fund.TestExogenousFund.java
/** * Test whether an instance of {@link ExogenousFund} behaves as * expected. This test operates as follows:<br><br> * /*from ww w. j a va 2s . co m*/ * {@code (a)} * An {@link ExogenousFund} entity is created. The deposits of this * {@link Fund} are held by an instance of {@link ExogenousDepositHolder};<br> * {@code (b)} * The {@link ExogenousFund} is provided with {@code floor(t)} units of * exogenous cash at simulation time {@code t}, for all {@code t};<br> * {@code (c)} * One {@link DepositingHousehold} {@link Agent} opens an investment * account with the {@link ExogenousFund}. {@code 1.0} is invested at * the time the account is created;<br> * {@code (d)} * Several simulation cycles elapse. After each cycle has elapsed, it is * asserted that the exogenous fund income has been transferred to the * {@link DepositingHousehold}. */ @Test public void testExogenousFundWithExogenousIncome() { final ExogenousDepositHolder depositHolder = new ExogenousDepositHolder(); final ExogenousFund fund = new ExogenousFund(depositHolder, new FundFactory() { @Override public Fund create(final DepositHolder depositHolder) { return new TrivialFund(depositHolder, new ClearingHouse()); } }, new FunctionParameter(new UnivariateFunction() { @Override public double value(double t) { return Math.floor(t); } }, "Parameter")); final DepositingHousehold household = new DepositingHousehold(depositHolder); household.debit(1.0); try { fund.invest(household, 1.); } catch (final InsufficientFundsException neverThrows) { Assert.fail(); } advanceTimeByOneCycle(); advanceTimeByOneCycle(); // Fires at T = 0 advanceTimeByOneCycle(); // Fires at T = 1 // At this point the exogenous fund has been debitted with 1.0 // from exogenous sources. This income should have been passed to // clients. Assert.assertEquals(household.getTotalAssets(), 2.0, 1.e-12); advanceTimeByOneCycle(); // Fires at T = 2 Assert.assertEquals(household.getTotalAssets(), 4.0, 1.e-12); advanceTimeByOneCycle(); // Fires at T = 3 Assert.assertEquals(household.getTotalAssets(), 7.0, 1.e-12); Assert.assertTrue(household.getTotalLiabilities() == 0.); Assert.assertEquals(fund.getTotalLiabilities(), 7.0, 1.e-12); Assert.assertEquals(fund.getTotalAssets(), 7.0, 1.e-12); Assert.assertEquals(fund.getBalance(household), 7.0, 1.e-12); }
From source file:de.unijena.bioinf.FragmentationTreeConstruction.computation.recalibration.AbstractRecalibrationStrategy.java
@Override public UnivariateFunction recalibrate(MutableSpectrum<Peak> spectrum, Spectrum<Peak> referenceSpectrum) { spectrum = new SimpleMutableSpectrum(spectrum); final SimpleMutableSpectrum ref = new SimpleMutableSpectrum(referenceSpectrum); preprocess(spectrum, ref);/* w ww . j a va 2 s . com*/ final double[][] values = MzRecalibration.maxIntervalStabbing(spectrum, referenceSpectrum, new UnivariateFunction() { @Override public double value(double x) { return epsilon.absoluteFor(x); } }, threshold); if (values[0].length < minNumberOfPeaks) return new Identity(); final UnivariateFunction recalibration = MzRecalibration.getMedianLinearRecalibration(values[0], values[1]); MzRecalibration.recalibrate(spectrum, recalibration); return recalibration; }
From source file:gedi.util.math.stat.distributions.GeneralizedExtremeValueDistribution.java
public static GeneralizedExtremeValueDistribution fitDataByMoments(double[] data, int start, int end) { int n = end - start; if (n < 2) throw new RuntimeException("Too few data!"); Arrays.sort(data, start, end); // compute moments final double[] b = new double[3]; for (int j = 1; j <= n; j++) { int index = j - 1 - start; b[0] += data[index];/* w w w. j av a 2 s . co m*/ b[1] += data[index] * ((j - 1.0) / (n - 1.0)); b[2] += data[index] * ((j - 1.0) / (n - 1.0)) * ((j - 2.0) / (n - 2.0)); } b[0] /= n; b[1] /= n; b[2] /= n; // solve parameters UnivariateFunction f = new UnivariateFunction() { public double value(double k) { return (3 * b[2] - b[0]) / (2 * b[1] - b[0]) - (1 - Math.pow(3, -k)) / (1 - Math.pow(2, -k)); } }; double k; if (Math.signum(f.value(-0.5)) != Math.signum(f.value(0.5))) k = UnivariateSolverUtils.solve(f, -0.5, 0.5); else { double c = (2 * b[1] - b[0]) / (3 * b[2] - b[0]) - Math.log(2) / Math.log(3); k = 7.859 * c + 2.9554 * c * c; } double g = Gamma.gamma(1 + k); double alpha = ((2 * b[1] - b[0]) * k) / (g * (1 - Math.pow(2, -k))); double xi = b[0] + alpha * (g - 1) / k; double location = xi; double scale = alpha; double shape = -k; return new GeneralizedExtremeValueDistribution(location, scale, shape); }
From source file:edu.rice.cs.bioinfo.programs.phylonet.algos.network.NetworkLikelihoodFromGTT.java
protected double findOptimalBranchLength(final Network<Object> speciesNetwork, final Map<String, List<String>> species2alleles, final List distinctTrees, final List gtCorrespondence, final Set<String> singleAlleleSpecies) { boolean continueRounds = true; // keep trying to improve network for (NetNode<Object> node : speciesNetwork.dfs()) { for (NetNode<Object> parent : node.getParents()) { node.setParentDistance(parent, 1.0); if (node.isNetworkNode()) { node.setParentProbability(parent, 0.5); }/*from w w w. ja v a 2s.c o m*/ } } Set<NetNode> node2ignoreForBL = findEdgeHavingNoBL(speciesNetwork, singleAlleleSpecies); double initalProb = computeProbabilityForCached(speciesNetwork, distinctTrees, species2alleles, gtCorrespondence); if (_printDetails) System.out.println(speciesNetwork.toString() + " : " + initalProb); final Container<Double> lnGtProbOfSpeciesNetwork = new Container<Double>(initalProb); // records the GTProb of the network at all times int roundIndex = 0; for (; roundIndex < _maxRounds && continueRounds; roundIndex++) { /* * Prepare a random ordering of network edge examinations each of which attempts to change a branch length or hybrid prob to improve the GTProb score. */ double lnGtProbLastRound = lnGtProbOfSpeciesNetwork.getContents(); List<Proc> assigmentActions = new ArrayList<Proc>(); // store adjustment commands here. Will execute them one by one later. for (final NetNode<Object> parent : edu.rice.cs.bioinfo.programs.phylonet.structs.network.util.Networks .postTraversal(speciesNetwork)) { for (final NetNode<Object> child : parent.getChildren()) { if (node2ignoreForBL.contains(child)) { continue; } assigmentActions.add(new Proc() { public void execute() { UnivariateFunction functionToOptimize = new UnivariateFunction() { public double value(double suggestedBranchLength) { double incumbentBranchLength = child.getParentDistance(parent); child.setParentDistance(parent, suggestedBranchLength); double lnProb = updateProbabilityForCached(speciesNetwork, distinctTrees, gtCorrespondence, child, parent); //System.out.println(speciesNetwork + ": " + lnProb); if (lnProb > lnGtProbOfSpeciesNetwork.getContents()) // did improve, keep change { lnGtProbOfSpeciesNetwork.setContents(lnProb); } else // didn't improve, roll back change { child.setParentDistance(parent, incumbentBranchLength); } return lnProb; } }; BrentOptimizer optimizer = new BrentOptimizer(_Brent1, _Brent2); // very small numbers so we control when brent stops, not brent. try { optimizer.optimize(_maxTryPerBranch, functionToOptimize, GoalType.MAXIMIZE, Double.MIN_VALUE, _maxBranchLength); } catch (TooManyEvaluationsException e) // _maxAssigmentAttemptsPerBranchParam exceeded { } updateProbabilityForCached(speciesNetwork, distinctTrees, gtCorrespondence, child, parent); if (_printDetails) System.out.println( speciesNetwork.toString() + " : " + lnGtProbOfSpeciesNetwork.getContents()); } }); } } for (final NetNode<Object> child : speciesNetwork.getNetworkNodes()) // find every hybrid node { Iterator<NetNode<Object>> hybridParents = child.getParents().iterator(); final NetNode hybridParent1 = hybridParents.next(); final NetNode hybridParent2 = hybridParents.next(); assigmentActions.add(new Proc() { public void execute() { UnivariateFunction functionToOptimize = new UnivariateFunction() { public double value(double suggestedProb) { double incumbentHybridProbParent1 = child.getParentProbability(hybridParent1); child.setParentProbability(hybridParent1, suggestedProb); child.setParentProbability(hybridParent2, 1.0 - suggestedProb); double lnProb = updateProbabilityForCached(speciesNetwork, distinctTrees, gtCorrespondence, child, null); //System.out.println(speciesNetwork + ": " + lnProb); if (lnProb > lnGtProbOfSpeciesNetwork.getContents()) // change improved GTProb, keep it { lnGtProbOfSpeciesNetwork.setContents(lnProb); } else // change did not improve, roll back { child.setParentProbability(hybridParent1, incumbentHybridProbParent1); child.setParentProbability(hybridParent2, 1.0 - incumbentHybridProbParent1); } return lnProb; } }; BrentOptimizer optimizer = new BrentOptimizer(_Brent1, _Brent2); // very small numbers so we control when brent stops, not brent. try { if (child.getName().equals("Y")) optimizer.optimize(_maxTryPerBranch, functionToOptimize, GoalType.MAXIMIZE, 0.6, 0.8); else optimizer.optimize(_maxTryPerBranch, functionToOptimize, GoalType.MAXIMIZE, 0, 1.0); } catch (TooManyEvaluationsException e) // _maxAssigmentAttemptsPerBranchParam exceeded { } updateProbabilityForCached(speciesNetwork, distinctTrees, gtCorrespondence, child, null); if (_printDetails) System.out.println( speciesNetwork.toString() + " : " + lnGtProbOfSpeciesNetwork.getContents()); } }); } // add hybrid probs to hybrid edges Collections.shuffle(assigmentActions); for (Proc assigment : assigmentActions) // for each change attempt, perform attempt { assigment.execute(); } if (_printDetails) { System.out.println("Round end ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"); System.out .println(speciesNetwork.toString() + "\n" + lnGtProbOfSpeciesNetwork.getContents() + "\n"); } if (((double) lnGtProbOfSpeciesNetwork.getContents()) == lnGtProbLastRound) // if no improvement was made wrt to last around, stop trying to find a better assignment { continueRounds = false; } else if (lnGtProbOfSpeciesNetwork.getContents() > lnGtProbLastRound) // improvement was made, ensure it is large enough wrt to improvement threshold to continue searching { double improvementPercentage = Math.pow(Math.E, (lnGtProbOfSpeciesNetwork.getContents() - lnGtProbLastRound)) - 1.0; // how much did we improve over last round if (improvementPercentage < _improvementThreshold) // improved, but not enough to keep searching { continueRounds = false; } } else { throw new IllegalStateException("Should never have decreased prob."); } } //System.out.println("\n" + lnGtProbOfSpeciesNetwork.getContents() + ": " + speciesNetwork); return lnGtProbOfSpeciesNetwork.getContents(); }
From source file:edu.uchc.octane.GaussianFitAstigmatism.java
/** * calculate the Z coordinate based on astigmatism *///from ww w. j a va 2 s .co m void calculateZ() { if (calibration_ == null) { z_ = 0; return; } UnivariateFunction func = new UnivariateFunction() { @Override public double value(double z) { double sigmax = FastMath.sqrt(pvp_.getPoint()[3] / 2); double sigmay = FastMath.sqrt(pvp_.getPoint()[4] / 2); double vx = calibration_[0] + (z - calibration_[1]) * (z - calibration_[1]) * calibration_[2] - sigmax; double vy = calibration_[3] + (z - calibration_[4]) * (z - calibration_[4]) * calibration_[5] - sigmay; return vx * vx + vy * vy; } }; BrentOptimizer optimizer = new BrentOptimizer(1e-4, 1e-4); UnivariatePointValuePair upvp = optimizer.optimize(new UnivariateObjectiveFunction(func), GoalType.MINIMIZE, MaxEval.unlimited(), new SearchInterval(2 * z0min_ - z0max_, 2 * z0max_ - z0min_)); if (upvp.getValue() > errTol_) { throw (new ConvergenceException()); } z_ = upvp.getPoint(); }
From source file:edu.utexas.cs.tactex.subscriptionspredictors.LWRCustOldAppache.java
/** * LWR prediction/*from w w w . ja va 2s .c o m*/ * * @param X * @param Y * @param x0 * @param tau * @return */ public Double LWRPredict(ArrayRealVector X, ArrayRealVector Y, double x0, final double tau) { ArrayRealVector X0 = new ArrayRealVector(X.getDimension(), x0); ArrayRealVector delta = X.subtract(X0); ArrayRealVector sqDists = delta.ebeMultiply(delta); UnivariateFunction expTau = new UnivariateFunction() { @Override public double value(double arg0) { //log.info(" cp univariate tau " + tau); return Math.pow(Math.E, -arg0 / (2 * tau)); } }; ArrayRealVector W = sqDists.map(expTau); double Xt_W_X = X.dotProduct(W.ebeMultiply(X)); if (Xt_W_X == 0.0) { log.error(" cp LWR cannot predict - 0 denominator returning NULL"); log.error("Xcv is " + X.toString()); log.error("Ycv is " + Y.toString()); log.error("x0 is " + x0); return null; // <==== NOTE: a caller must be prepared for it } double theta = (1.0 / Xt_W_X) * X.ebeMultiply(W).dotProduct(Y); return theta * x0; }