1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.apache.commons.math.distribution;
18
19 import java.io.Serializable;
20
21 import org.apache.commons.math.MathException;
22 import org.apache.commons.math.util.MathUtils;
23
24 /***
25 * The default implementation of {@link HypergeometricDistribution}.
26 *
27 * @version $Revision: 1.16 $ $Date: 2004/11/07 03:32:48 $
28 */
29 public class HypergeometricDistributionImpl extends AbstractIntegerDistribution
30 implements HypergeometricDistribution, Serializable
31 {
32
33 /*** Serializable version identifier */
34 static final long serialVersionUID = -436928820673516179L;
35
36 /*** The number of successes in the population. */
37 private int numberOfSuccesses;
38
39 /*** The population size. */
40 private int populationSize;
41
42 /*** The sample size. */
43 private int sampleSize;
44
45 /***
46 * Construct a new hypergeometric distribution with the given the population
47 * size, the number of successes in the population, and the sample size.
48 * @param populationSize the population size.
49 * @param numberOfSuccesses number of successes in the population.
50 * @param sampleSize the sample size.
51 */
52 public HypergeometricDistributionImpl(int populationSize,
53 int numberOfSuccesses, int sampleSize) {
54 super();
55 if (numberOfSuccesses > populationSize) {
56 throw new IllegalArgumentException(
57 "number of successes must be less than or equal to population size");
58 }
59 if (sampleSize > populationSize) {
60 throw new IllegalArgumentException(
61 "sample size must be less than or equal to population size");
62 }
63 setPopulationSize(populationSize);
64 setSampleSize(sampleSize);
65 setNumberOfSuccesses(numberOfSuccesses);
66 }
67
68 /***
69 * For this disbution, X, this method returns P(X ≤ x).
70 * @param x the value at which the PDF is evaluated.
71 * @return PDF for this distribution.
72 * @throws MathException if the cumulative probability can not be
73 * computed due to convergence or other numerical errors.
74 */
75 public double cumulativeProbability(int x) throws MathException{
76 double ret;
77
78 int n = getPopulationSize();
79 int m = getNumberOfSuccesses();
80 int k = getSampleSize();
81
82 int[] domain = getDomain(n, m, k);
83 if (x < domain[0]) {
84 ret = 0.0;
85 } else if(x >= domain[1]) {
86 ret = 1.0;
87 } else {
88 ret = 0.0;
89 for (int i = domain[0]; i <= x; ++i){
90 ret += probability(n, m, k, i);
91 }
92 }
93
94 return ret;
95 }
96
97 /***
98 * Return the domain for the given hypergeometric distribution parameters.
99 * @param n the population size.
100 * @param m number of successes in the population.
101 * @param k the sample size.
102 * @return a two element array containing the lower and upper bounds of the
103 * hypergeometric distribution.
104 */
105 private int[] getDomain(int n, int m, int k){
106 return new int[]{
107 getLowerDomain(n, m, k),
108 getUpperDomain(m, k)
109 };
110 }
111
112 /***
113 * Access the domain value lower bound, based on <code>p</code>, used to
114 * bracket a PDF root.
115 *
116 * @param p the desired probability for the critical value
117 * @return domain value lower bound, i.e.
118 * P(X < <i>lower bound</i>) < <code>p</code>
119 */
120 protected int getDomainLowerBound(double p) {
121 return getLowerDomain(getPopulationSize(), getNumberOfSuccesses(),
122 getSampleSize());
123 }
124
125 /***
126 * Access the domain value upper bound, based on <code>p</code>, used to
127 * bracket a PDF root.
128 *
129 * @param p the desired probability for the critical value
130 * @return domain value upper bound, i.e.
131 * P(X < <i>upper bound</i>) > <code>p</code>
132 */
133 protected int getDomainUpperBound(double p) {
134 return getUpperDomain(getSampleSize(), getNumberOfSuccesses());
135 }
136
137 /***
138 * Return the lowest domain value for the given hypergeometric distribution
139 * parameters.
140 * @param n the population size.
141 * @param m number of successes in the population.
142 * @param k the sample size.
143 * @return the lowest domain value of the hypergeometric distribution.
144 */
145 private int getLowerDomain(int n, int m, int k) {
146 return Math.max(0, m - (n - k));
147 }
148
149 /***
150 * Access the number of successes.
151 * @return the number of successes.
152 */
153 public int getNumberOfSuccesses() {
154 return numberOfSuccesses;
155 }
156
157 /***
158 * Access the population size.
159 * @return the population size.
160 */
161 public int getPopulationSize() {
162 return populationSize;
163 }
164
165 /***
166 * Access the sample size.
167 * @return the sample size.
168 */
169 public int getSampleSize() {
170 return sampleSize;
171 }
172
173 /***
174 * Return the highest domain value for the given hypergeometric distribution
175 * parameters.
176 * @param m number of successes in the population.
177 * @param k the sample size.
178 * @return the highest domain value of the hypergeometric distribution.
179 */
180 private int getUpperDomain(int m, int k){
181 return Math.min(k, m);
182 }
183
184 /***
185 * For this disbution, X, this method returns P(X = x).
186 *
187 * @param x the value at which the PMF is evaluated.
188 * @return PMF for this distribution.
189 */
190 public double probability(int x) {
191 double ret;
192
193 int n = getPopulationSize();
194 int m = getNumberOfSuccesses();
195 int k = getSampleSize();
196
197 int[] domain = getDomain(n, m, k);
198 if(x < domain[0] || x > domain[1]){
199 ret = 0.0;
200 } else {
201 ret = probability(n, m, k, x);
202 }
203
204 return ret;
205 }
206
207 /***
208 * For the disbution, X, defined by the given hypergeometric distribution
209 * parameters, this method returns P(X = x).
210 *
211 * @param n the population size.
212 * @param m number of successes in the population.
213 * @param k the sample size.
214 * @param x the value at which the PMF is evaluated.
215 * @return PMF for the distribution.
216 */
217 private double probability(int n, int m, int k, int x) {
218 return Math.exp(MathUtils.binomialCoefficientLog(m, x) +
219 MathUtils.binomialCoefficientLog(n - m, k - x) -
220 MathUtils.binomialCoefficientLog(n, k));
221 }
222
223 /***
224 * Modify the number of successes.
225 * @param num the new number of successes.
226 * @throws IllegalArgumentException if <code>num</code> is negative.
227 */
228 public void setNumberOfSuccesses(int num) {
229 if(num < 0){
230 throw new IllegalArgumentException(
231 "number of successes must be non-negative.");
232 }
233 numberOfSuccesses = num;
234 }
235
236 /***
237 * Modify the population size.
238 * @param size the new population size.
239 * @throws IllegalArgumentException if <code>size</code> is not positive.
240 */
241 public void setPopulationSize(int size) {
242 if(size <= 0){
243 throw new IllegalArgumentException(
244 "population size must be positive.");
245 }
246 populationSize = size;
247 }
248
249 /***
250 * Modify the sample size.
251 * @param size the new sample size.
252 * @throws IllegalArgumentException if <code>size</code> is negative.
253 */
254 public void setSampleSize(int size) {
255 if (size < 0) {
256 throw new IllegalArgumentException(
257 "sample size must be non-negative.");
258 }
259 sampleSize = size;
260 }
261 }