001 /*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements. See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License. You may obtain a copy of the License at
008 *
009 * http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017 package org.apache.commons.math.distribution;
018
019 import java.io.Serializable;
020
021
022 /**
023 * The default implementation of {@link ChiSquaredDistribution}
024 *
025 * @version $Id: ChiSquaredDistributionImpl.java 1178295 2011-10-03 04:36:27Z psteitz $
026 */
027 public class ChiSquaredDistributionImpl
028 extends AbstractContinuousDistribution
029 implements ChiSquaredDistribution, Serializable {
030 /**
031 * Default inverse cumulative probability accuracy
032 * @since 2.1
033 */
034 public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9;
035 /** Serializable version identifier */
036 private static final long serialVersionUID = -8352658048349159782L;
037 /** Internal Gamma distribution. */
038 private final GammaDistribution gamma;
039 /** Inverse cumulative probability accuracy */
040 private final double solverAbsoluteAccuracy;
041
042 /**
043 * Create a Chi-Squared distribution with the given degrees of freedom.
044 *
045 * @param degreesOfFreedom Degrees of freedom.
046 */
047 public ChiSquaredDistributionImpl(double degreesOfFreedom) {
048 this(degreesOfFreedom, DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
049 }
050
051 /**
052 * Create a Chi-Squared distribution with the given degrees of freedom and
053 * inverse cumulative probability accuracy.
054 *
055 * @param degreesOfFreedom Degrees of freedom.
056 * @param inverseCumAccuracy the maximum absolute error in inverse
057 * cumulative probability estimates (defaults to
058 * {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}).
059 * @since 2.1
060 */
061 public ChiSquaredDistributionImpl(double degreesOfFreedom,
062 double inverseCumAccuracy) {
063 gamma = new GammaDistributionImpl(degreesOfFreedom / 2, 2);
064 solverAbsoluteAccuracy = inverseCumAccuracy;
065 }
066
067 /**
068 * {@inheritDoc}
069 */
070 public double getDegreesOfFreedom() {
071 return gamma.getAlpha() * 2.0;
072 }
073
074 /**
075 * {@inheritDoc}
076 */
077 @Override
078 public double density(double x) {
079 return gamma.density(x);
080 }
081
082 /**
083 * For this distribution, {@code X}, this method returns {@code P(X < x)}.
084 *
085 * @param x the value at which the CDF is evaluated.
086 * @return CDF for this distribution.
087 */
088 public double cumulativeProbability(double x) {
089 return gamma.cumulativeProbability(x);
090 }
091
092 /**
093 * For this distribution, X, this method returns the critical point
094 * {@code x}, such that {@code P(X < x) = p}.
095 * It will return 0 when p = 0 and {@code Double.POSITIVE_INFINITY}
096 * when p = 1.
097 *
098 * @param p Desired probability.
099 * @return {@code x}, such that {@code P(X < x) = p}.
100 * @throws org.apache.commons.math.exception.OutOfRangeException if
101 * {@code p} is not a valid probability.
102 */
103 @Override
104 public double inverseCumulativeProbability(final double p) {
105 if (p == 0) {
106 return 0d;
107 }
108 if (p == 1) {
109 return Double.POSITIVE_INFINITY;
110 }
111 return super.inverseCumulativeProbability(p);
112 }
113
114 /**
115 * Access the domain value lower bound, based on {@code p}, used to
116 * bracket a CDF root. This method is used by
117 * {@link #inverseCumulativeProbability(double)} to find critical values.
118 *
119 * @param p the desired probability for the critical value
120 * @return domain value lower bound, i.e. {@code P(X < 'lower bound') < p}.
121 */
122 @Override
123 protected double getDomainLowerBound(double p) {
124 return Double.MIN_VALUE * gamma.getBeta();
125 }
126
127 /**
128 * Access the domain value upper bound, based on {@code p}, used to
129 * bracket a CDF root. This method is used by
130 * {@link #inverseCumulativeProbability(double)} to find critical values.
131 *
132 * @param p Desired probability for the critical value.
133 * @return domain value upper bound, i.e. {@code P(X < 'upper bound') > p}.
134 */
135 @Override
136 protected double getDomainUpperBound(double p) {
137 // NOTE: chi squared is skewed to the left
138 // NOTE: therefore, P(X < μ) > .5
139
140 double ret;
141
142 if (p < .5) {
143 // use mean
144 ret = getDegreesOfFreedom();
145 } else {
146 // use max
147 ret = Double.MAX_VALUE;
148 }
149
150 return ret;
151 }
152
153 /**
154 * Access the initial domain value, based on {@code p}, used to
155 * bracket a CDF root. This method is used by
156 * {@link #inverseCumulativeProbability(double)} to find critical values.
157 *
158 * @param p Desired probability for the critical value.
159 * @return the initial domain value.
160 */
161 @Override
162 protected double getInitialDomain(double p) {
163 // NOTE: chi squared is skewed to the left
164 // NOTE: therefore, P(X < μ) > 0.5
165
166 double ret;
167
168 if (p < 0.5) {
169 // use 1/2 mean
170 ret = getDegreesOfFreedom() * 0.5;
171 } else {
172 // use mean
173 ret = getDegreesOfFreedom();
174 }
175
176 return ret;
177 }
178
179 /**
180 * Return the absolute accuracy setting of the solver used to estimate
181 * inverse cumulative probabilities.
182 *
183 * @return the solver absolute accuracy.
184 * @since 2.1
185 */
186 @Override
187 protected double getSolverAbsoluteAccuracy() {
188 return solverAbsoluteAccuracy;
189 }
190
191 /**
192 * {@inheritDoc}
193 *
194 * The lower bound of the support is always 0 no matter the
195 * degrees of freedom.
196 *
197 * @return lower bound of the support (always 0)
198 */
199 @Override
200 public double getSupportLowerBound() {
201 return 0;
202 }
203
204 /**
205 * {@inheritDoc}
206 *
207 * The upper bound of the support is always positive infinity no matter the
208 * degrees of freedom.
209 *
210 * @return upper bound of the support (always Double.POSITIVE_INFINITY)
211 */
212 @Override
213 public double getSupportUpperBound() {
214 return Double.POSITIVE_INFINITY;
215 }
216
217 /**
218 * {@inheritDoc}
219 *
220 * For <code>k</code> degrees of freedom, the mean is
221 * <code>k</code>
222 *
223 * @return {@inheritDoc}
224 */
225 @Override
226 protected double calculateNumericalMean() {
227 return getDegreesOfFreedom();
228 }
229
230 /**
231 * {@inheritDoc}
232 *
233 * For <code>k</code> degrees of freedom, the variance is
234 * <code>2 * k</code>
235 *
236 * @return {@inheritDoc}
237 */
238 @Override
239 protected double calculateNumericalVariance() {
240 return 2*getDegreesOfFreedom();
241 }
242
243 /**
244 * {@inheritDoc}
245 */
246 @Override
247 public boolean isSupportLowerBoundInclusive() {
248 return true;
249 }
250
251 /**
252 * {@inheritDoc}
253 */
254 @Override
255 public boolean isSupportUpperBoundInclusive() {
256 return false;
257 }
258 }