source: src/main/java/agents/anac/y2019/harddealer/math3/distribution/FDistribution.java

Last change on this file was 204, checked in by Katsuhide Fujita, 5 years ago

Fixed errors of ANAC2019 agents
  • Property svn:executable set to *
File size: 11.5 KB
Line 
1/*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17
18package agents.anac.y2019.harddealer.math3.distribution;
19
20import agents.anac.y2019.harddealer.math3.exception.NotStrictlyPositiveException;
21import agents.anac.y2019.harddealer.math3.exception.util.LocalizedFormats;
22import agents.anac.y2019.harddealer.math3.random.RandomGenerator;
23import agents.anac.y2019.harddealer.math3.random.Well19937c;
24import agents.anac.y2019.harddealer.math3.special.Beta;
25import agents.anac.y2019.harddealer.math3.util.FastMath;
26
27/**
28 * Implementation of the F-distribution.
29 *
30 * @see <a href="http://en.wikipedia.org/wiki/F-distribution">F-distribution (Wikipedia)</a>
31 * @see <a href="http://mathworld.wolfram.com/F-Distribution.html">F-distribution (MathWorld)</a>
32 */
33public class FDistribution extends AbstractRealDistribution {
34 /**
35 * Default inverse cumulative probability accuracy.
36 * @since 2.1
37 */
38 public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9;
39 /** Serializable version identifier. */
40 private static final long serialVersionUID = -8516354193418641566L;
41 /** The numerator degrees of freedom. */
42 private final double numeratorDegreesOfFreedom;
43 /** The numerator degrees of freedom. */
44 private final double denominatorDegreesOfFreedom;
45 /** Inverse cumulative probability accuracy. */
46 private final double solverAbsoluteAccuracy;
47 /** Cached numerical variance */
48 private double numericalVariance = Double.NaN;
49 /** Whether or not the numerical variance has been calculated */
50 private boolean numericalVarianceIsCalculated = false;
51
52 /**
53 * Creates an F distribution using the given degrees of freedom.
54 * <p>
55 * <b>Note:</b> this constructor will implicitly create an instance of
56 * {@link Well19937c} as random generator to be used for sampling only (see
57 * {@link #sample()} and {@link #sample(int)}). In case no sampling is
58 * needed for the created distribution, it is advised to pass {@code null}
59 * as random generator via the appropriate constructors to avoid the
60 * additional initialisation overhead.
61 *
62 * @param numeratorDegreesOfFreedom Numerator degrees of freedom.
63 * @param denominatorDegreesOfFreedom Denominator degrees of freedom.
64 * @throws NotStrictlyPositiveException if
65 * {@code numeratorDegreesOfFreedom <= 0} or
66 * {@code denominatorDegreesOfFreedom <= 0}.
67 */
68 public FDistribution(double numeratorDegreesOfFreedom,
69 double denominatorDegreesOfFreedom)
70 throws NotStrictlyPositiveException {
71 this(numeratorDegreesOfFreedom, denominatorDegreesOfFreedom,
72 DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
73 }
74
75 /**
76 * Creates an F distribution using the given degrees of freedom
77 * and inverse cumulative probability accuracy.
78 * <p>
79 * <b>Note:</b> this constructor will implicitly create an instance of
80 * {@link Well19937c} as random generator to be used for sampling only (see
81 * {@link #sample()} and {@link #sample(int)}). In case no sampling is
82 * needed for the created distribution, it is advised to pass {@code null}
83 * as random generator via the appropriate constructors to avoid the
84 * additional initialisation overhead.
85 *
86 * @param numeratorDegreesOfFreedom Numerator degrees of freedom.
87 * @param denominatorDegreesOfFreedom Denominator degrees of freedom.
88 * @param inverseCumAccuracy the maximum absolute error in inverse
89 * cumulative probability estimates.
90 * @throws NotStrictlyPositiveException if
91 * {@code numeratorDegreesOfFreedom <= 0} or
92 * {@code denominatorDegreesOfFreedom <= 0}.
93 * @since 2.1
94 */
95 public FDistribution(double numeratorDegreesOfFreedom,
96 double denominatorDegreesOfFreedom,
97 double inverseCumAccuracy)
98 throws NotStrictlyPositiveException {
99 this(new Well19937c(), numeratorDegreesOfFreedom,
100 denominatorDegreesOfFreedom, inverseCumAccuracy);
101 }
102
103 /**
104 * Creates an F distribution.
105 *
106 * @param rng Random number generator.
107 * @param numeratorDegreesOfFreedom Numerator degrees of freedom.
108 * @param denominatorDegreesOfFreedom Denominator degrees of freedom.
109 * @throws NotStrictlyPositiveException if {@code numeratorDegreesOfFreedom <= 0} or
110 * {@code denominatorDegreesOfFreedom <= 0}.
111 * @since 3.3
112 */
113 public FDistribution(RandomGenerator rng,
114 double numeratorDegreesOfFreedom,
115 double denominatorDegreesOfFreedom)
116 throws NotStrictlyPositiveException {
117 this(rng, numeratorDegreesOfFreedom, denominatorDegreesOfFreedom, DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
118 }
119
120 /**
121 * Creates an F distribution.
122 *
123 * @param rng Random number generator.
124 * @param numeratorDegreesOfFreedom Numerator degrees of freedom.
125 * @param denominatorDegreesOfFreedom Denominator degrees of freedom.
126 * @param inverseCumAccuracy the maximum absolute error in inverse
127 * cumulative probability estimates.
128 * @throws NotStrictlyPositiveException if {@code numeratorDegreesOfFreedom <= 0} or
129 * {@code denominatorDegreesOfFreedom <= 0}.
130 * @since 3.1
131 */
132 public FDistribution(RandomGenerator rng,
133 double numeratorDegreesOfFreedom,
134 double denominatorDegreesOfFreedom,
135 double inverseCumAccuracy)
136 throws NotStrictlyPositiveException {
137 super(rng);
138
139 if (numeratorDegreesOfFreedom <= 0) {
140 throw new NotStrictlyPositiveException(LocalizedFormats.DEGREES_OF_FREEDOM,
141 numeratorDegreesOfFreedom);
142 }
143 if (denominatorDegreesOfFreedom <= 0) {
144 throw new NotStrictlyPositiveException(LocalizedFormats.DEGREES_OF_FREEDOM,
145 denominatorDegreesOfFreedom);
146 }
147 this.numeratorDegreesOfFreedom = numeratorDegreesOfFreedom;
148 this.denominatorDegreesOfFreedom = denominatorDegreesOfFreedom;
149 solverAbsoluteAccuracy = inverseCumAccuracy;
150 }
151
152 /**
153 * {@inheritDoc}
154 *
155 * @since 2.1
156 */
157 public double density(double x) {
158 return FastMath.exp(logDensity(x));
159 }
160
161 /** {@inheritDoc} **/
162 @Override
163 public double logDensity(double x) {
164 final double nhalf = numeratorDegreesOfFreedom / 2;
165 final double mhalf = denominatorDegreesOfFreedom / 2;
166 final double logx = FastMath.log(x);
167 final double logn = FastMath.log(numeratorDegreesOfFreedom);
168 final double logm = FastMath.log(denominatorDegreesOfFreedom);
169 final double lognxm = FastMath.log(numeratorDegreesOfFreedom * x +
170 denominatorDegreesOfFreedom);
171 return nhalf * logn + nhalf * logx - logx +
172 mhalf * logm - nhalf * lognxm - mhalf * lognxm -
173 Beta.logBeta(nhalf, mhalf);
174 }
175
176 /**
177 * {@inheritDoc}
178 *
179 * The implementation of this method is based on
180 * <ul>
181 * <li>
182 * <a href="http://mathworld.wolfram.com/F-Distribution.html">
183 * F-Distribution</a>, equation (4).
184 * </li>
185 * </ul>
186 */
187 public double cumulativeProbability(double x) {
188 double ret;
189 if (x <= 0) {
190 ret = 0;
191 } else {
192 double n = numeratorDegreesOfFreedom;
193 double m = denominatorDegreesOfFreedom;
194
195 ret = Beta.regularizedBeta((n * x) / (m + n * x),
196 0.5 * n,
197 0.5 * m);
198 }
199 return ret;
200 }
201
202 /**
203 * Access the numerator degrees of freedom.
204 *
205 * @return the numerator degrees of freedom.
206 */
207 public double getNumeratorDegreesOfFreedom() {
208 return numeratorDegreesOfFreedom;
209 }
210
211 /**
212 * Access the denominator degrees of freedom.
213 *
214 * @return the denominator degrees of freedom.
215 */
216 public double getDenominatorDegreesOfFreedom() {
217 return denominatorDegreesOfFreedom;
218 }
219
220 /** {@inheritDoc} */
221 @Override
222 protected double getSolverAbsoluteAccuracy() {
223 return solverAbsoluteAccuracy;
224 }
225
226 /**
227 * {@inheritDoc}
228 *
229 * For denominator degrees of freedom parameter {@code b}, the mean is
230 * <ul>
231 * <li>if {@code b > 2} then {@code b / (b - 2)},</li>
232 * <li>else undefined ({@code Double.NaN}).
233 * </ul>
234 */
235 public double getNumericalMean() {
236 final double denominatorDF = getDenominatorDegreesOfFreedom();
237
238 if (denominatorDF > 2) {
239 return denominatorDF / (denominatorDF - 2);
240 }
241
242 return Double.NaN;
243 }
244
245 /**
246 * {@inheritDoc}
247 *
248 * For numerator degrees of freedom parameter {@code a} and denominator
249 * degrees of freedom parameter {@code b}, the variance is
250 * <ul>
251 * <li>
252 * if {@code b > 4} then
253 * {@code [2 * b^2 * (a + b - 2)] / [a * (b - 2)^2 * (b - 4)]},
254 * </li>
255 * <li>else undefined ({@code Double.NaN}).
256 * </ul>
257 */
258 public double getNumericalVariance() {
259 if (!numericalVarianceIsCalculated) {
260 numericalVariance = calculateNumericalVariance();
261 numericalVarianceIsCalculated = true;
262 }
263 return numericalVariance;
264 }
265
266 /**
267 * used by {@link #getNumericalVariance()}
268 *
269 * @return the variance of this distribution
270 */
271 protected double calculateNumericalVariance() {
272 final double denominatorDF = getDenominatorDegreesOfFreedom();
273
274 if (denominatorDF > 4) {
275 final double numeratorDF = getNumeratorDegreesOfFreedom();
276 final double denomDFMinusTwo = denominatorDF - 2;
277
278 return ( 2 * (denominatorDF * denominatorDF) * (numeratorDF + denominatorDF - 2) ) /
279 ( (numeratorDF * (denomDFMinusTwo * denomDFMinusTwo) * (denominatorDF - 4)) );
280 }
281
282 return Double.NaN;
283 }
284
285 /**
286 * {@inheritDoc}
287 *
288 * The lower bound of the support is always 0 no matter the parameters.
289 *
290 * @return lower bound of the support (always 0)
291 */
292 public double getSupportLowerBound() {
293 return 0;
294 }
295
296 /**
297 * {@inheritDoc}
298 *
299 * The upper bound of the support is always positive infinity
300 * no matter the parameters.
301 *
302 * @return upper bound of the support (always Double.POSITIVE_INFINITY)
303 */
304 public double getSupportUpperBound() {
305 return Double.POSITIVE_INFINITY;
306 }
307
308 /** {@inheritDoc} */
309 public boolean isSupportLowerBoundInclusive() {
310 return false;
311 }
312
313 /** {@inheritDoc} */
314 public boolean isSupportUpperBoundInclusive() {
315 return false;
316 }
317
318 /**
319 * {@inheritDoc}
320 *
321 * The support of this distribution is connected.
322 *
323 * @return {@code true}
324 */
325 public boolean isSupportConnected() {
326 return true;
327 }
328}
Note: See TracBrowser for help on using the repository browser.