1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18 package org.apache.commons.math.distribution;
19
20 import java.io.Serializable;
21
22 import org.apache.commons.math.MathRuntimeException;
23 import org.apache.commons.math.util.MathUtils;
24
25
26
27
28
29
30 public class HypergeometricDistributionImpl extends AbstractIntegerDistribution
31 implements HypergeometricDistribution, Serializable
32 {
33
34
35 private static final long serialVersionUID = -436928820673516179L;
36
37
38 private int numberOfSuccesses;
39
40
41 private int populationSize;
42
43
44 private int sampleSize;
45
46
47
48
49
50
51
52
53 public HypergeometricDistributionImpl(int populationSize,
54 int numberOfSuccesses, int sampleSize) {
55 super();
56 if (numberOfSuccesses > populationSize) {
57 throw MathRuntimeException.createIllegalArgumentException(
58 "number of successes ({0}) must be less than or equal to population size ({1})",
59 numberOfSuccesses, populationSize);
60 }
61 if (sampleSize > populationSize) {
62 throw MathRuntimeException.createIllegalArgumentException(
63 "sample size ({0}) must be less than or equal to population size ({1})",
64 sampleSize, populationSize);
65 }
66 setPopulationSize(populationSize);
67 setSampleSize(sampleSize);
68 setNumberOfSuccesses(numberOfSuccesses);
69 }
70
71
72
73
74
75
76 @Override
77 public double cumulativeProbability(int x) {
78 double ret;
79
80 int n = getPopulationSize();
81 int m = getNumberOfSuccesses();
82 int k = getSampleSize();
83
84 int[] domain = getDomain(n, m, k);
85 if (x < domain[0]) {
86 ret = 0.0;
87 } else if(x >= domain[1]) {
88 ret = 1.0;
89 } else {
90 ret = innerCumulativeProbability(domain[0], x, 1, n, m, k);
91 }
92
93 return ret;
94 }
95
96
97
98
99
100
101
102
103
104 private int[] getDomain(int n, int m, int k){
105 return new int[]{
106 getLowerDomain(n, m, k),
107 getUpperDomain(m, k)
108 };
109 }
110
111
112
113
114
115
116
117
118
119 @Override
120 protected int getDomainLowerBound(double p) {
121 return getLowerDomain(getPopulationSize(), getNumberOfSuccesses(),
122 getSampleSize());
123 }
124
125
126
127
128
129
130
131
132
133 @Override
134 protected int getDomainUpperBound(double p) {
135 return getUpperDomain(getSampleSize(), getNumberOfSuccesses());
136 }
137
138
139
140
141
142
143
144
145
146 private int getLowerDomain(int n, int m, int k) {
147 return Math.max(0, m - (n - k));
148 }
149
150
151
152
153
154 public int getNumberOfSuccesses() {
155 return numberOfSuccesses;
156 }
157
158
159
160
161
162 public int getPopulationSize() {
163 return populationSize;
164 }
165
166
167
168
169
170 public int getSampleSize() {
171 return sampleSize;
172 }
173
174
175
176
177
178
179
180
181 private int getUpperDomain(int m, int k){
182 return Math.min(k, m);
183 }
184
185
186
187
188
189
190
191 public double probability(int x) {
192 double ret;
193
194 int n = getPopulationSize();
195 int m = getNumberOfSuccesses();
196 int k = getSampleSize();
197
198 int[] domain = getDomain(n, m, k);
199 if(x < domain[0] || x > domain[1]){
200 ret = 0.0;
201 } else {
202 ret = probability(n, m, k, x);
203 }
204
205 return ret;
206 }
207
208
209
210
211
212
213
214
215
216
217
218 private double probability(int n, int m, int k, int x) {
219 return Math.exp(MathUtils.binomialCoefficientLog(m, x) +
220 MathUtils.binomialCoefficientLog(n - m, k - x) -
221 MathUtils.binomialCoefficientLog(n, k));
222 }
223
224
225
226
227
228
229 public void setNumberOfSuccesses(int num) {
230 if(num < 0){
231 throw MathRuntimeException.createIllegalArgumentException(
232 "number of successes must be non-negative ({0})",
233 num);
234 }
235 numberOfSuccesses = num;
236 }
237
238
239
240
241
242
243 public void setPopulationSize(int size) {
244 if(size <= 0){
245 throw MathRuntimeException.createIllegalArgumentException(
246 "population size must be positive ({0})",
247 size);
248 }
249 populationSize = size;
250 }
251
252
253
254
255
256
257 public void setSampleSize(int size) {
258 if (size < 0) {
259 throw MathRuntimeException.createIllegalArgumentException(
260 "sample size must be positive ({0})",
261 size);
262 }
263 sampleSize = size;
264 }
265
266
267
268
269
270
271
272 public double upperCumulativeProbability(int x) {
273 double ret;
274
275 int n = getPopulationSize();
276 int m = getNumberOfSuccesses();
277 int k = getSampleSize();
278
279 int[] domain = getDomain(n, m, k);
280 if (x < domain[0]) {
281 ret = 1.0;
282 } else if(x > domain[1]) {
283 ret = 0.0;
284 } else {
285 ret = innerCumulativeProbability(domain[1], x, -1, n, m, k);
286 }
287
288 return ret;
289 }
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304 private double innerCumulativeProbability(
305 int x0, int x1, int dx, int n, int m, int k)
306 {
307 double ret = probability(n, m, k, x0);
308 while (x0 != x1) {
309 x0 += dx;
310 ret += probability(n, m, k, x0);
311 }
312 return ret;
313 }
314 }