001: /*
002: * <copyright>
003: *
004: * Copyright 1997-2004 BBNT Solutions, LLC
005: * under sponsorship of the Defense Advanced Research Projects
006: * Agency (DARPA).
007: *
008: * You can redistribute this software and/or modify it under the
009: * terms of the Cougaar Open Source License as published on the
010: * Cougaar Open Source Website (www.cougaar.org).
011: *
012: * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
013: * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
014: * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
015: * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
016: * OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
017: * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
018: * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
019: * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
020: * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
021: * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
022: * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
023: *
024: * </copyright>
025: */
026: package org.cougaar.util;
027:
028: import java.util.Iterator;
029: import java.util.SortedMap;
030: import java.util.TreeMap;
031:
032: public class Random extends java.util.Random {
033: // Memoize for repeated use with the same mean rate (xm)
034: private double oldXm = -1.0; // Flag for whether xm has changed since last call,
035: private double g, sq, logXm;
036:
037: /**
038: * Returns an integer value that is a random deviate drawn from a
039: * Poisson distribution of mean xm.
040: **/
041: public int nextPoisson(double xm) {
042: double em, t;
043: if (xm < 25.0) { // Use direct method (25 is the empirically
044: if (xm != oldXm) { // determined performance crossover).
045: oldXm = xm;
046: g = Math.exp(-xm); // If Xta is new, compute the exponential.
047: }
048: em = -1.0;
049: t = 1.0;
050: while (true) {
051: em += 1.0; // Instead of adding exponential deviates it is equivalent to
052: t *= nextDouble(); // multiply uniform deviates. Then we never actually have to
053: if (t <= g)
054: return (int) em; // take the log, merely compare to the pre-computed exponential.
055: }
056: } else { // Use rejection method.
057: if (xm != oldXm) { // If xm has changed since the last call, then precompute
058: oldXm = xm; // some functions which occur below.
059: sq = Math.sqrt(2.0 * xm);
060: logXm = Math.log(xm);
061: g = xm * logXm - // The function gammaLn is
062: MoreMath.gammaLn(xm + 1.0); // the natural log of the gamma function, as given in 6.2.
063: }
064: double y;
065: while (true) {
066: y = Math.tan(Math.PI * nextDouble()); // y is a deviate from a Lorentzian comparison function.
067: em = sq * y + xm; // em is y, shifted and scaled.
068: if (em < 0.0)
069: continue; // Reject if in regime of zero probability.
070: em = Math.floor(em); // The trick for integer-valued distributions.
071: t = (0.9 // the ratio of the desired distribution to
072: * (1.0 + y * y) // the comparison function;
073: * Math.exp(em * logXm // we accept or reject by comparing it
074: - MoreMath.gammaLn(em + 1.0) - g)); // to another uniform deviate. The factor 0.9 is chosen
075: if (nextDouble() > t)
076: continue; // so that t never exceeds 1.
077: break;
078: }
079: }
080: return (int) em;
081: }
082:
083: /* Test harness. Generate a set of samples with specified mean */
084: /*
085: public static void main(String[] args) {
086: double xm = 1.0;
087: double nTrials = 100;
088: if (args.length > 0) xm = Double.parseDouble(args[0]);
089: if (args.length > 1) nTrials = Double.parseDouble(args[1]);
090: SortedMap bins = new TreeMap();
091: Random random = new Random();
092: class Bin {
093: int count;
094: };
095: double sum = 0.0;
096: for (int i = 0; i < nTrials; i++) {
097: int q = random.nextPoisson(xm);
098: Integer k = new Integer(q);
099: Bin bin = (Bin) bins.get(k);
100: if (bin == null) {
101: bin = new Bin();
102: bins.put(k, bin);
103: }
104: bin.count++;
105: sum += q;
106: }
107: for (Iterator i = bins.keySet().iterator(); i.hasNext(); ) {
108: Integer k = (Integer) i.next();
109: Bin bin = (Bin) bins.get(k);
110: System.out.println(k + "\t" + bin.count);
111: }
112: }
113: */
114: }
|