1 package org.lcsim.math.distribution;
2
3 import java.util.Random;
4
5 /**
6 * Generate random numbers according to an empirical distribution provided by
7 * the user as an array of positive real numbers.
8 * Based on RandomGeneral from CLHEP.
9 *
10 * @author Norman A. Graf
11 *
12 * @version $Id:
13 */
14 public class EmpiricalDistribution
15 {
16 private Random _rand;
17 protected double[] _cdf; // cumulative distribution function
18 protected int _interpolationType;
19
20 /**
21 * Create an EmpiricalDistribution from an input array of positive values.
22 * The array value i should be thought of as the central value of bin i in
23 * an equally binned distribution. The values will be normalized internally.
24 * This constructor will return doubles from the usual interval between
25 * 0 and 1.
26 * @param pdf an array of positive values representing the empirical
27 * distribution from which random numbers should be drawn.
28 */
29 public EmpiricalDistribution(double[] pdf)
30 {
31 _rand = new Random();
32 setState(pdf);
33 }
34
35 //TODO implement method which allows random number range to be
36 //set by user beforehand
37 // public EmpiricalDistribution(double[] pdf, double lowEdge, double interval)
38 // {
39 // _rand = new Random();
40 // setState(pdf);
41 // }
42
43 /**
44 * Allows the random number generator seed to be set
45 * @param seed
46 */
47 public void setSeed(long seed)
48 {
49 _rand.setSeed(seed);
50 }
51
52 /**
53 * Returns a random number from the distribution.
54 * @return the next double value drawn from this distribution
55 */
56 public double nextDouble()
57 {
58 double rand = _rand.nextDouble();
59 if (this._cdf == null)
60 {
61 return rand; // Non-existing pdf
62 }
63
64 int nBins = _cdf.length - 1;
65 int nbelow = 0; // largest k such that I[k] is known to be <= rand
66 int nabove = nBins; // largest k such that I[k] is known to be > rand
67
68 while (nabove > nbelow + 1)
69 {
70 int middle = (nabove + nbelow + 1) >> 1; // div 2
71 if (rand >= _cdf[middle])
72 {
73 nbelow = middle;
74 } else
75 {
76 nabove = middle;
77 }
78 }
79 // nabove is always nbelow+1 and they straddle rand:
80 double binMeasure = _cdf[nabove] - _cdf[nbelow];
81 if (binMeasure == 0.0)
82 {
83 // rand lies right in a bin of measure 0.
84 // Return the center of the range of that bin.
85 // (Any value between k/N and (k+1)/N is equally good,
86 // in this rare case.)
87 return (nbelow + 0.5) / nBins;
88 }
89 double binFraction = (rand - _cdf[nbelow]) / binMeasure;
90 return (nbelow + binFraction) / nBins;
91 }
92
93 /**
94 * Returns the probability distribution function.
95 */
96 public double pdf(int k)
97 {
98 if (k < 0 || k >= _cdf.length - 1)
99 {
100 return 0.0;
101 }
102 return _cdf[k] - _cdf[k - 1];
103 }
104
105 /**
106 * Creates and normalizes the cumulative probability distribution
107 */
108 public void setState(double[] pdf)
109 {
110 if (pdf == null || pdf.length == 0)
111 {
112 this._cdf = null;
113 //throw new IllegalArgumentException("Non-existing pdf");
114 return;
115 }
116 // compute cumulative distribution function (cdf) from probability
117 // distribution function (pdf)
118 int nBins = pdf.length;
119 this._cdf = new double[nBins + 1];
120
121 _cdf[0] = 0;
122 for (int ptn = 0; ptn < nBins; ++ptn)
123 {
124 double prob = pdf[ptn];
125 if (prob < 0.0)
126 {
127 throw new IllegalArgumentException("Negative probability");
128 }
129 _cdf[ptn + 1] = _cdf[ptn] + prob;
130 }
131 if (_cdf[nBins] <= 0.0)
132 {
133 throw new IllegalArgumentException("At least one probability must be > 0.0");
134 }
135 //normalize to 1.
136 for (int ptn = 0; ptn < nBins + 1; ++ptn)
137 {
138 _cdf[ptn] /= _cdf[nBins];
139 }
140
141 }
142
143 /**
144 * Returns a String representation of the receiver.
145 */
146 public String toString()
147 {
148 return this.getClass().getName();
149 }
150 }