1 /*
2 Copyright � 1999 CERN - European Organization for Nuclear Research.
3 Permission to use, copy, modify, distribute and sell this software and its documentation for any purpose
4 is hereby granted without fee, provided that the above copyright notice appear in all copies and
5 that both that copyright notice and this permission notice appear in supporting documentation.
6 CERN makes no representations about the suitability of this software for any purpose.
7 It is provided "as is" without expressed or implied warranty.
8 */
9 package cern.jet.random;
10
11 import cern.jet.random.engine.RandomEngine;
12 /**
13 * Logarithmic distribution.
14 * <p>
15 * Valid parameter ranges: <tt>0 < p < 1</tt>.
16 * <p>
17 * Instance methods operate on a user supplied uniform random number generator; they are unsynchronized.
18 * <dt>
19 * Static methods operate on a default uniform random number generator; they are synchronized.
20 * <p>
21 * <b>Implementation:</b>
22 * <dt>
23 * Method: Inversion/Transformation.
24 * <dt>
25 * This is a port of <tt>lsk.c</tt> from the <A HREF="http://www.cis.tu-graz.ac.at/stat/stadl/random.html">C-RAND / WIN-RAND</A> library.
26 * C-RAND's implementation, in turn, is based upon
27 * <p>
28 * A.W. Kemp (1981): Efficient generation of logarithmically distributed pseudo-random variables, Appl. Statist. 30, 249-253.
29 *
30 * @author wolfgang.hoschek@cern.ch
31 * @version 1.0, 09/24/99
32 */
33 public class Logarithmic extends AbstractContinousDistribution {
34 protected double my_p;
35
36 // cached vars for method nextDouble(a) (for performance only)
37 private double t,h,a_prev = -1.0;
38
39 // The uniform random number generated shared by all <b>static</b> methods.
40 protected static Logarithmic shared = new Logarithmic(0.5,makeDefaultGenerator());
41 /**
42 * Constructs a Logarithmic distribution.
43 */
Logarithmic(double p, RandomEngine randomGenerator)44 public Logarithmic(double p, RandomEngine randomGenerator) {
45 setRandomGenerator(randomGenerator);
46 setState(p);
47 }
48 /**
49 * Returns a random number from the distribution.
50 */
nextDouble()51 public double nextDouble() {
52 return nextDouble(this.my_p);
53 }
54 /**
55 * Returns a random number from the distribution; bypasses the internal state.
56 */
nextDouble(double a)57 public double nextDouble(double a) {
58 /******************************************************************
59 * *
60 * Logarithmic Distribution - Inversion/Transformation *
61 * *
62 ******************************************************************
63 * *
64 * The algorithm combines Inversion and Transformation. *
65 * It is based on the following fact: A random variable X from *
66 * the Logarithmic distribution has the property that X for fixed *
67 * Y=y is Geometric distributed with P(X=x|Y=y)=(1-y)*y^(x-1) (*) *
68 * where Y has distribution function F(y)=ln(1-y)/ln(1-p). *
69 * So first random numbers y are generated by simple Inversion, *
70 * then k=(long int) (1+ln(u)/ln(y)) is a Geometric random number *
71 * and because of (*) a Logarithmic one. *
72 * To speed up the algorithm squeezes are used as well as the *
73 * fact, that many of the random numbers are 1 or 2 (depending on *
74 * special circumstances). *
75 * On an IBM/PC 486 optimal performance is achieved, if for p<.97 *
76 * simple inversion is used and otherwise the transformation. *
77 * On an IBM/PC 286 inversion should be restricted to p<.90. *
78 * *
79 ******************************************************************
80 * *
81 * FUNCTION: - lsk samples a random number from the *
82 * Logarithmic distribution with *
83 * parameter 0 < p < 1 . *
84 * REFERENCE: - A.W. Kemp (1981): Efficient generation of *
85 * logarithmically distributed pseudo-random *
86 * variables, Appl. Statist. 30, 249-253. *
87 * SUBPROGRAMS: - drand(seed) ... (0,1)-Uniform generator with *
88 * unsigned long integer *seed. *
89 * *
90 ******************************************************************/
91 double u,v,p,q;
92 int k;
93
94 if (a != a_prev) { // Set-up
95 a_prev = a;
96 if (a<0.97) t = -a / Math.log(1.0 - a);
97 else h=Math.log(1.0 - a);
98 }
99
100 u=randomGenerator.raw();
101 if (a<0.97) { // Inversion/Chop-down
102 k = 1;
103 p = t;
104 while (u > p) {
105 //System.out.println("u="+u+", p="+p);
106 u -= p;
107 k++;
108 p *= a * (k-1.0)/(double)k;
109 }
110 return k;
111 }
112
113 if (u > a) return 1; // Transformation
114 u=randomGenerator.raw();
115 v = u;
116 q = 1.0 - Math.exp(v * h);
117 if ( u <= q * q) {
118 k = (int) (1 + Math.log(u) / Math.log(q));
119 return k;
120 }
121 if (u > q) return 1;
122 return 2;
123 }
124 /**
125 * Sets the distribution parameter.
126 */
setState(double p)127 public void setState(double p) {
128 this.my_p = p;
129 }
130 /**
131 * Returns a random number from the distribution.
132 */
staticNextDouble(double p)133 public static double staticNextDouble(double p) {
134 synchronized (shared) {
135 return shared.nextDouble(p);
136 }
137 }
138 /**
139 * Returns a String representation of the receiver.
140 */
toString()141 public String toString() {
142 return this.getClass().getName()+"("+my_p+")";
143 }
144 /**
145 * Sets the uniform random number generated shared by all <b>static</b> methods.
146 * @param randomGenerator the new uniform random number generator to be shared.
147 */
xstaticSetRandomGenerator(RandomEngine randomGenerator)148 private static void xstaticSetRandomGenerator(RandomEngine randomGenerator) {
149 synchronized (shared) {
150 shared.setRandomGenerator(randomGenerator);
151 }
152 }
153 }
154