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
018 package org.apache.commons.math.distribution;
019
020 import java.io.Serializable;
021
022 import org.apache.commons.math.exception.NotStrictlyPositiveException;
023 import org.apache.commons.math.exception.OutOfRangeException;
024 import org.apache.commons.math.exception.util.LocalizedFormats;
025 import org.apache.commons.math.util.FastMath;
026
027 /**
028 * Default implementation of
029 * {@link org.apache.commons.math.distribution.CauchyDistribution}.
030 *
031 * @since 1.1
032 * @version $Id: CauchyDistributionImpl.java 1131229 2011-06-03 20:49:25Z luc $
033 */
034 public class CauchyDistributionImpl extends AbstractContinuousDistribution
035 implements CauchyDistribution, Serializable {
036 /**
037 * Default inverse cumulative probability accuracy.
038 * @since 2.1
039 */
040 public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9;
041 /** Serializable version identifier */
042 private static final long serialVersionUID = 8589540077390120676L;
043 /** The median of this distribution. */
044 private final double median;
045 /** The scale of this distribution. */
046 private final double scale;
047 /** Inverse cumulative probability accuracy */
048 private final double solverAbsoluteAccuracy;
049
050 /**
051 * Creates cauchy distribution with the medain equal to zero and scale
052 * equal to one.
053 */
054 public CauchyDistributionImpl() {
055 this(0, 1);
056 }
057
058 /**
059 * Create a cauchy distribution using the given median and scale.
060 *
061 * @param median Median for this distribution.
062 * @param scale Scale parameter for this distribution.
063 */
064 public CauchyDistributionImpl(double median, double scale) {
065 this(median, scale, DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
066 }
067
068 /**
069 * Create a cauchy distribution using the given median and scale.
070 *
071 * @param median Median for this distribution.
072 * @param scale Scale parameter for this distribution.
073 * @param inverseCumAccuracy Maximum absolute error in inverse
074 * cumulative probability estimates
075 * (defaults to {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}).
076 * @throws NotStrictlyPositiveException if {@code s <= 0}.
077 * @since 2.1
078 */
079 public CauchyDistributionImpl(double median, double scale,
080 double inverseCumAccuracy) {
081 if (scale <= 0) {
082 throw new NotStrictlyPositiveException(LocalizedFormats.SCALE, scale);
083 }
084 this.scale = scale;
085 this.median = median;
086 solverAbsoluteAccuracy = inverseCumAccuracy;
087 }
088
089 /**
090 * For this distribution, {@code X}, this method returns {@code P(X < x)}.
091 *
092 * @param x Value at which the CDF is evaluated.
093 * @return CDF evaluated at {@code x}.
094 */
095 public double cumulativeProbability(double x) {
096 return 0.5 + (FastMath.atan((x - median) / scale) / FastMath.PI);
097 }
098
099 /**
100 * {@inheritDoc}
101 */
102 public double getMedian() {
103 return median;
104 }
105
106 /**
107 * {@inheritDoc}
108 */
109 public double getScale() {
110 return scale;
111 }
112
113 /**
114 * {@inheritDoc}
115 */
116 @Override
117 public double density(double x) {
118 final double dev = x - median;
119 return (1 / FastMath.PI) * (scale / (dev * dev + scale * scale));
120 }
121
122 /**
123 * For this distribution, {@code X}, this method returns the critical
124 * point {@code x}, such that {@code P(X < x) = p}.
125 * It will return {@code Double.NEGATIVE_INFINITY} when p = 0 and
126 * {@code Double.POSITIVE_INFINITY} when p = 1.
127 *
128 * @param p Desired probability.
129 * @return {@code x}, such that {@code P(X < x) = p}.
130 * @throws OutOfRangeException if {@code p} is not a valid probability.
131 */
132 @Override
133 public double inverseCumulativeProbability(double p) {
134 double ret;
135 if (p < 0 || p > 1) {
136 throw new OutOfRangeException(p, 0, 1);
137 } else if (p == 0) {
138 ret = Double.NEGATIVE_INFINITY;
139 } else if (p == 1) {
140 ret = Double.POSITIVE_INFINITY;
141 } else {
142 ret = median + scale * FastMath.tan(FastMath.PI * (p - .5));
143 }
144 return ret;
145 }
146
147 /**
148 * Access the domain value lower bound, based on {@code p}, used to
149 * bracket a CDF root. This method is used by
150 * {@link #inverseCumulativeProbability(double)} to find critical values.
151 *
152 * @param p Desired probability for the critical value.
153 * @return domain value lower bound, i.e. {@code P(X < 'lower bound') < p}.
154 */
155 @Override
156 protected double getDomainLowerBound(double p) {
157 double ret;
158
159 if (p < 0.5) {
160 ret = -Double.MAX_VALUE;
161 } else {
162 ret = median;
163 }
164
165 return ret;
166 }
167
168 /**
169 * Access the domain value upper bound, based on <code>p</code>, used to
170 * bracket a CDF root. This method is used by
171 * {@link #inverseCumulativeProbability(double)} to find critical values.
172 *
173 * @param p Desired probability for the critical value.
174 * @return domain value lower bound, i.e. {@code P(X < 'upper bound') > p}.
175 */
176 @Override
177 protected double getDomainUpperBound(double p) {
178 double ret;
179
180 if (p < 0.5) {
181 ret = median;
182 } else {
183 ret = Double.MAX_VALUE;
184 }
185
186 return ret;
187 }
188
189 /**
190 * Access the initial domain value, based on {@code p}, used to
191 * bracket a CDF root. This method is used by
192 * {@link #inverseCumulativeProbability(double)} to find critical values.
193 *
194 * @param p Desired probability for the critical value.
195 * @return the initial domain value.
196 */
197 @Override
198 protected double getInitialDomain(double p) {
199 double ret;
200
201 if (p < 0.5) {
202 ret = median - scale;
203 } else if (p > 0.5) {
204 ret = median + scale;
205 } else {
206 ret = median;
207 }
208
209 return ret;
210 }
211
212 /**
213 * Return the absolute accuracy setting of the solver used to estimate
214 * inverse cumulative probabilities.
215 *
216 * @return the solver absolute accuracy
217 * @since 2.1
218 */
219 @Override
220 protected double getSolverAbsoluteAccuracy() {
221 return solverAbsoluteAccuracy;
222 }
223
224 /**
225 * {@inheritDoc}
226 *
227 * The lower bound of the support is always negative infinity no matter
228 * the parameters.
229 *
230 * @return lower bound of the support (always Double.NEGATIVE_INFINITY)
231 */
232 @Override
233 public double getSupportLowerBound() {
234 return Double.NEGATIVE_INFINITY;
235 }
236
237 /**
238 * {@inheritDoc}
239 *
240 * The upper bound of the support is always positive infinity no matter
241 * the parameters.
242 *
243 * @return upper bound of the support (always Double.POSITIVE_INFINITY)
244 */
245 @Override
246 public double getSupportUpperBound() {
247 return Double.POSITIVE_INFINITY;
248 }
249
250 /**
251 * {@inheritDoc}
252 *
253 * The mean is always undefined no matter the parameters.
254 *
255 * @return mean (always Double.NaN)
256 */
257 @Override
258 protected double calculateNumericalMean() {
259 return Double.NaN;
260 }
261
262 /**
263 * {@inheritDoc}
264 *
265 * The variance is always undefined no matter the parameters.
266 *
267 * @return variance (always Double.NaN)
268 */
269 @Override
270 protected double calculateNumericalVariance() {
271 return Double.NaN;
272 }
273
274 /**
275 * {@inheritDoc}
276 */
277 @Override
278 public boolean isSupportLowerBoundInclusive() {
279 return false;
280 }
281
282 /**
283 * {@inheritDoc}
284 */
285 @Override
286 public boolean isSupportUpperBoundInclusive() {
287 return false;
288 }
289 }