List of usage examples for org.apache.commons.math.util FastMath log
public static double log(final double x)
From source file:cz.paulrz.montecarlo.single.LogArrivedPointValuation.java
/** {@inheritedDoc} */ public final Double value(final Path path) { final double[] values = path.getValues(); return FastMath.log(values[path.getLength() - 1] / values[0]); }
From source file:cz.paulrz.montecarlo.single.ExpOrnsteinUhlenbeckProcess.java
/** {@inheritDoc} */ @Override/*from w ww. j a v a2 s . c om*/ public double drift(double t, double x) { return (theta * (mu - FastMath.log(x)) + halfsigmasquare) * x; }
From source file:net.malariagen.gatk.math.SaddlePointExpansion.java
/** * Compute the error of Stirling's series at the given value. * <p>//w ww . j a va2 s .c o m * References: * <ol> * <li>Eric W. Weisstein. "Stirling's Series." From MathWorld--A Wolfram Web * Resource. <a target="_blank" * href="http://mathworld.wolfram.com/StirlingsSeries.html"> * http://mathworld.wolfram.com/StirlingsSeries.html</a></li> * </ol> * </p> * * @param z the value. * @return the Striling's series error. */ static double getStirlingError(double z) { double ret; if (z < 15.0) { double z2 = 2.0 * z; if (FastMath.floor(z2) == z2) { ret = EXACT_STIRLING_ERRORS[(int) z2]; } else { ret = Gamma.logGamma(z + 1.0) - (z + 0.5) * FastMath.log(z) + z - HALF_LOG_2_PI; } } else { double z2 = z * z; ret = (0.083333333333333333333 - (0.00277777777777777777778 - (0.00079365079365079365079365 - (0.000595238095238095238095238 - 0.0008417508417508417508417508 / z2) / z2) / z2) / z2) / z; } return ret; }
From source file:net.malariagen.gatk.math.SaddlePointExpansion.java
/** * A part of the deviance portion of the saddle point approximation. * <p>/*ww w . j a v a2 s.c o m*/ * References: * <ol> * <li>Catherine Loader (2000). "Fast and Accurate Computation of Binomial * Probabilities.". <a target="_blank" * href="http://www.herine.net/stat/papers/dbinom.pdf"> * http://www.herine.net/stat/papers/dbinom.pdf</a></li> * </ol> * </p> * * @param x the x value. * @param mu the average. * @return a part of the deviance. */ static double getDeviancePart(double x, double mu) { double ret; if (FastMath.abs(x - mu) < 0.1 * (x + mu)) { double d = x - mu; double v = d / (x + mu); double s1 = v * d; double s = Double.NaN; double ej = 2.0 * x * v; v = v * v; int j = 1; while (s1 != s) { s = s1; ej *= v; s1 = s + ej / ((j * 2) + 1); ++j; } ret = s1; } else { ret = x * FastMath.log(x / mu) + mu - x; } return ret; }
From source file:net.malariagen.gatk.math.SaddlePointExpansion.java
/** * Compute the PMF for a binomial distribution using the saddle point * expansion.//from w w w . j av a 2s. c o m * * @param x the value at which the probability is evaluated. * @param n the number of trials. * @param p the probability of success. * @param q the probability of failure (1 - p). * @return log(p(x)). */ public static double logBinomialProbability(int x, int n, double p, double q) { double ret; if (x == 0) { if (p < 0.1) { ret = -getDeviancePart(n, n * q) - n * p; } else { ret = n * FastMath.log(q); } } else if (x == n) { if (q < 0.1) { ret = -getDeviancePart(n, n * p) - n * q; } else { ret = n * FastMath.log(p); } } else { ret = getStirlingError(n) - getStirlingError(x) - getStirlingError(n - x) - getDeviancePart(x, n * p) - getDeviancePart(n - x, n * q); double f = (MathUtils.TWO_PI * x * (n - x)) / n; ret = -0.5 * FastMath.log(f) + ret; } return ret; }
From source file:cz.paulrz.montecarlo.random.Sobol.java
@Override public double nextGaussian() { final double random; if (Double.isNaN(nextGaussian)) { // generate a new pair of gaussian numbers final double[] xs = nextPoint(); final double x = xs[0]; final double y = xs[1]; final double alpha = 2 * FastMath.PI * x; final double r = FastMath.sqrt(-2 * FastMath.log(y)); random = r * FastMath.cos(alpha); nextGaussian = r * FastMath.sin(alpha); } else {// w ww .j av a2s . c om // use the second element of the pair already generated random = nextGaussian; nextGaussian = Double.NaN; } return random; }
From source file:org.onebusaway.nyc.vehicle_tracking.impl.inference.likelihood.ScheduleLikelihood.java
/** * Manual truncation (when using the formal run distribution) * We shouldn't have to worry about normalization, yet. *//* ww w. ja va 2 s . c om*/ public static double getScheduleDevLogProb(final double x, StudentTDistribution schedDist) { final double pSched; final boolean isFormal = schedDist == schedDevFormalRunDist; if (!isFormal || (x <= POS_SCHED_DEV_CUTOFF && x >= NEG_SCHED_DEV_CUTOFF)) { pSched = (isFormal ? FastMath.log(pFormal) : FastMath.log1p(-pFormal)) + schedDist.getProbabilityFunction().logEvaluate(x); } else { pSched = Double.NEGATIVE_INFINITY; } return pSched; }
From source file:org.onebusaway.nyc.vehicle_tracking.impl.inference.MotionModelImpl.java
@Override public Multiset<Particle> move(Multiset<Particle> particles, double timestamp, double timeElapsed, Observation obs, boolean previouslyResampled) throws ParticleFilterException { final Multiset<Particle> results = HashMultiset.create(); boolean anySnapped = false; /*//www. j a v a2 s. com * TODO FIXME could keep a particle-set-info record... */ for (final Multiset.Entry<Particle> parent : particles.entrySet()) { final VehicleState state = parent.getElement().getData(); if (state.getBlockStateObservation() != null && state.getBlockStateObservation().isSnapped()) { anySnapped = true; break; } } final double ess = ParticleFilter.getEffectiveSampleSize(particles); final boolean generalBlockTransition = allowGeneralBlockTransition(obs, ess, previouslyResampled, anySnapped); double normOffset = Double.NEGATIVE_INFINITY; for (final Multiset.Entry<Particle> parent : particles.entrySet()) { final VehicleState parentState = parent.getElement().getData(); final BlockStateObservation parentBlockStateObs = parentState.getBlockStateObservation(); final Set<BlockStateObservation> transitions = Sets.newHashSet(); if (generalBlockTransition || allowParentBlockTransition(parentState, obs)) { /* * These are all the snapped and DSC/run blocks */ transitions.addAll(_blocksFromObservationService.determinePotentialBlockStatesForObservation(obs)); } else if (parentBlockStateObs != null) { /* * Only the snapped blocks. * We are also allowing changes to snapped in-progress states when * they were sampled well outside of allowed backward search distance. */ final double backwardDistance = Double.NEGATIVE_INFINITY; transitions.addAll(_blocksFromObservationService.advanceState(obs, parentState.getBlockState(), backwardDistance, Double.POSITIVE_INFINITY)); } BlockStateObservation newParentBlockStateObs; if (parentBlockStateObs != null) { newParentBlockStateObs = _blocksFromObservationService.getBlockStateObservationFromDist(obs, parentBlockStateObs.getBlockState().getBlockInstance(), parentBlockStateObs.getBlockState().getBlockLocation().getDistanceAlongBlock()); } else { newParentBlockStateObs = null; } transitions.add(newParentBlockStateObs); /* * We make a subtle distinction here, by allowing the previous state to * remain as a transition but unaltered. This helps in cases of * deadhead->in-progress transitions, since, later on, the block state we * create in the following will be treated differently, signaled by way of * pointer equality. */ final double vehicleHasNotMovedProb = MovedLikelihood.computeVehicleHasNotMovedProbability(obs); for (int i = 0; i < parent.getCount(); ++i) { final Particle sampledParticle = sampleTransitionParticle(parent, newParentBlockStateObs, obs, vehicleHasNotMovedProb, transitions); normOffset = LogMath.add(sampledParticle.getLogWeight(), normOffset); results.add(sampledParticle); } } /* * Normalize */ for (final Entry<Particle> p : results.entrySet()) { final double thisTotalWeight = p.getElement().getLogWeight() + FastMath.log(p.getCount()); double logNormed = thisTotalWeight - normOffset; if (logNormed > 0d) logNormed = 0d; p.getElement().setLogNormedWeight(logNormed); } return results; }
From source file:org.onebusaway.nyc.vehicle_tracking.impl.inference.ParticleFactoryImpl.java
@Override public Multiset<Particle> createParticles(double timestamp, Observation obs) throws ParticleFilterException { final Set<BlockStateObservation> potentialBlocks = _blocksFromObservationService .determinePotentialBlockStatesForObservation(obs); final Multiset<Particle> particles = HashMultiset.create(); double normOffset = Double.NEGATIVE_INFINITY; for (int i = 0; i < _initialNumberOfParticles; ++i) { final CategoricalDist<Particle> transitionProb = new CategoricalDist<Particle>(); for (final BlockStateObservation blockState : potentialBlocks) { final SensorModelResult transProb = new SensorModelResult("transition"); final double inMotionSample = threadLocalRng.get().nextDouble(); final boolean vehicleNotMoved = inMotionSample < 0.5; final MotionState motionState = _motionModel.updateMotionState(obs, vehicleNotMoved); BlockStateObservation sampledBlockState; if (blockState != null) { /*/*from w ww. j av a 2 s . c om*/ * Sample a distance along the block using the snapped observation * results as priors. */ if (blockState.isSnapped()) { sampledBlockState = blockState; } else { sampledBlockState = _blockStateSamplingStrategy .samplePriorScheduleState(blockState.getBlockState().getBlockInstance(), obs); } } else { sampledBlockState = null; } final JourneyState journeyState = _journeyStateTransitionModel.getJourneyState(sampledBlockState, null, obs, vehicleNotMoved); final VehicleState state = vehicleState(motionState, sampledBlockState, journeyState, obs); final Context context = new Context(null, state, obs); transProb.addResultAsAnd(_motionModel.getEdgeLikelihood().likelihood(context)); transProb.addResultAsAnd(_motionModel.getGpsLikelihood().likelihood(context)); transProb.addResultAsAnd(_motionModel.getSchedLikelihood().likelihood(context)); transProb.addResultAsAnd(_motionModel.dscLikelihood.likelihood(context)); transProb.addResultAsAnd(_motionModel.runLikelihood.likelihood(context)); transProb.addResultAsAnd(_motionModel.runTransitionLikelihood.likelihood(context)); transProb.addResultAsAnd(_motionModel.nullStateLikelihood.likelihood(context)); transProb.addResultAsAnd(_motionModel.nullLocationLikelihood.likelihood(context)); final Particle newParticle = new Particle(timestamp, null, 0.0, state); newParticle.setResult(transProb); transitionProb.logPut(transProb.getLogProbability(), newParticle); } final Particle newSample; if (transitionProb.canSample()) { newSample = transitionProb.sample(); newSample.setLogWeight(newSample.getResult().getLogProbability()); particles.add(newSample); } else { final double inMotionSample = ParticleFactoryImpl.getThreadLocalRng().get().nextDouble(); final boolean vehicleNotMoved = inMotionSample < 0.5; final MotionState motionState = _motionModel.updateMotionState(obs, vehicleNotMoved); final JourneyState journeyState = _journeyStateTransitionModel.getJourneyState(null, null, obs, vehicleNotMoved); final VehicleState nullState = new VehicleState(motionState, null, journeyState, null, obs); final Context context = new Context(null, nullState, obs); final SensorModelResult priorProb = new SensorModelResult("prior creation"); priorProb.addResultAsAnd(_motionModel.getEdgeLikelihood().likelihood(context)); priorProb.addResultAsAnd(_motionModel.getGpsLikelihood().likelihood(context)); priorProb.addResultAsAnd(_motionModel.getSchedLikelihood().likelihood(context)); priorProb.addResultAsAnd(_motionModel.dscLikelihood.likelihood(context)); priorProb.addResultAsAnd(_motionModel.runLikelihood.likelihood(context)); priorProb.addResultAsAnd(_motionModel.runTransitionLikelihood.likelihood(context)); priorProb.addResultAsAnd(_motionModel.nullStateLikelihood.likelihood(context)); priorProb.addResultAsAnd(_motionModel.nullLocationLikelihood.likelihood(context)); newSample = new Particle(timestamp, null, 0.0, nullState); newSample.setResult(priorProb); particles.add(newSample); newSample.setLogWeight(newSample.getResult().getLogProbability()); } normOffset = LogMath.add(newSample.getLogWeight(), normOffset); } /* * Normalize */ for (final Entry<Particle> p : particles.entrySet()) { p.getElement() .setLogNormedWeight(p.getElement().getLogWeight() + FastMath.log(p.getCount()) - normOffset); } return particles; }
From source file:org.onebusaway.nyc.vehicle_tracking.impl.particlefilter.CategoricalDistTest.java
@Test public void testSampleA() throws ParticleFilterException { final CategoricalDist<String> cdf = new CategoricalDist<String>(); cdf.logPut(FastMath.log(0.3 * 1e-5), "c"); cdf.logPut(FastMath.log(0.3), "c"); cdf.logPut(FastMath.log(0.2), "c"); cdf.logPut(FastMath.log(0.01), "a"); cdf.logPut(FastMath.log(0.001), "a"); cdf.logPut(FastMath.log(0.2 * 1e-7), "b"); final Counter<String> counter = new Counter<String>(); final int iterations = 1000; for (int i = 0; i < iterations; i++) counter.increment(cdf.sample()); final double a = counter.getCount("a") / (double) iterations; final double b = counter.getCount("b") / (double) iterations; final double c = counter.getCount("c") / (double) iterations; final double cummulativeProb = cdf.getCummulativeProb(); assertEquals(a, cdf.density("a") / cummulativeProb, .05); assertEquals(b, cdf.density("b") / cummulativeProb, .05); assertEquals(c, cdf.density("c") / cummulativeProb, .05); cdf.logPut(FastMath.log(0.001), "d"); final Multiset<String> res = cdf.sample(iterations); assertEquals(res.count("a") / (double) iterations, cdf.density("a") / cummulativeProb, .05); assertEquals(res.count("b") / (double) iterations, cdf.density("b") / cummulativeProb, .05); assertEquals(res.count("c") / (double) iterations, cdf.density("c") / cummulativeProb, .05); assertEquals(res.count("d") / (double) iterations, cdf.density("d") / cummulativeProb, .05); final Multiset<Particle> testSet = HashMultiset.create(); final Particle p1 = new Particle(0, null, cdf.density("a")); testSet.add(p1);/*from w w w .ja v a 2 s. c o m*/ final Particle p2 = new Particle(0, null, cdf.density("b")); testSet.add(p2); final Particle p3 = new Particle(0, null, cdf.density("c")); testSet.add(p3); final Particle p4 = new Particle(0, null, cdf.density("d")); testSet.add(p4); final Multiset<Particle> resSet = ParticleFilter.lowVarianceSampler(testSet, iterations); final Counter<Particle> counter2 = new Counter<Particle>(); for (final Particle p : resSet) counter2.increment(p); final double a2 = counter2.getCount(p1) / (double) iterations; final double b2 = counter2.getCount(p2) / (double) iterations; final double c2 = counter2.getCount(p3) / (double) iterations; final double d2 = counter2.getCount(p4) / (double) iterations; assertEquals(a2, cdf.density("a") / cummulativeProb, .05); assertEquals(b2, cdf.density("b") / cummulativeProb, .05); assertEquals(c2, cdf.density("c") / cummulativeProb, .05); assertEquals(d2, cdf.density("d") / cummulativeProb, .05); }