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