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 */
017package org.apache.commons.math3.distribution;
018
019import org.apache.commons.math3.exception.NotStrictlyPositiveException;
020import org.apache.commons.math3.exception.OutOfRangeException;
021import org.apache.commons.math3.exception.util.LocalizedFormats;
022import org.apache.commons.math3.random.RandomGenerator;
023import org.apache.commons.math3.random.Well19937c;
024import org.apache.commons.math3.util.FastMath;
025
026/**
027 * Implementation of the Cauchy distribution.
028 *
029 * @see <a href="http://en.wikipedia.org/wiki/Cauchy_distribution">Cauchy distribution (Wikipedia)</a>
030 * @see <a href="http://mathworld.wolfram.com/CauchyDistribution.html">Cauchy Distribution (MathWorld)</a>
031 * @since 1.1 (changed to concrete class in 3.0)
032 */
033public class CauchyDistribution extends AbstractRealDistribution {
034    /**
035     * Default inverse cumulative probability accuracy.
036     * @since 2.1
037     */
038    public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9;
039    /** Serializable version identifier */
040    private static final long serialVersionUID = 8589540077390120676L;
041    /** The median of this distribution. */
042    private final double median;
043    /** The scale of this distribution. */
044    private final double scale;
045    /** Inverse cumulative probability accuracy */
046    private final double solverAbsoluteAccuracy;
047
048    /**
049     * Creates a Cauchy distribution with the median equal to zero and scale
050     * equal to one.
051     */
052    public CauchyDistribution() {
053        this(0, 1);
054    }
055
056    /**
057     * Creates a Cauchy distribution using the given median and scale.
058     * <p>
059     * <b>Note:</b> this constructor will implicitly create an instance of
060     * {@link Well19937c} as random generator to be used for sampling only (see
061     * {@link #sample()} and {@link #sample(int)}). In case no sampling is
062     * needed for the created distribution, it is advised to pass {@code null}
063     * as random generator via the appropriate constructors to avoid the
064     * additional initialisation overhead.
065     *
066     * @param median Median for this distribution.
067     * @param scale Scale parameter for this distribution.
068     */
069    public CauchyDistribution(double median, double scale) {
070        this(median, scale, DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
071    }
072
073    /**
074     * Creates a Cauchy distribution using the given median and scale.
075     * <p>
076     * <b>Note:</b> this constructor will implicitly create an instance of
077     * {@link Well19937c} as random generator to be used for sampling only (see
078     * {@link #sample()} and {@link #sample(int)}). In case no sampling is
079     * needed for the created distribution, it is advised to pass {@code null}
080     * as random generator via the appropriate constructors to avoid the
081     * additional initialisation overhead.
082     *
083     * @param median Median for this distribution.
084     * @param scale Scale parameter for this distribution.
085     * @param inverseCumAccuracy Maximum absolute error in inverse
086     * cumulative probability estimates
087     * (defaults to {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}).
088     * @throws NotStrictlyPositiveException if {@code scale <= 0}.
089     * @since 2.1
090     */
091    public CauchyDistribution(double median, double scale,
092                              double inverseCumAccuracy) {
093        this(new Well19937c(), median, scale, inverseCumAccuracy);
094    }
095
096    /**
097     * Creates a Cauchy distribution.
098     *
099     * @param rng Random number generator.
100     * @param median Median for this distribution.
101     * @param scale Scale parameter for this distribution.
102     * @throws NotStrictlyPositiveException if {@code scale <= 0}.
103     * @since 3.3
104     */
105    public CauchyDistribution(RandomGenerator rng, double median, double scale) {
106        this(rng, median, scale, DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
107    }
108
109    /**
110     * Creates a Cauchy distribution.
111     *
112     * @param rng Random number generator.
113     * @param median Median for this distribution.
114     * @param scale Scale parameter for this distribution.
115     * @param inverseCumAccuracy Maximum absolute error in inverse
116     * cumulative probability estimates
117     * (defaults to {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}).
118     * @throws NotStrictlyPositiveException if {@code scale <= 0}.
119     * @since 3.1
120     */
121    public CauchyDistribution(RandomGenerator rng,
122                              double median,
123                              double scale,
124                              double inverseCumAccuracy) {
125        super(rng);
126        if (scale <= 0) {
127            throw new NotStrictlyPositiveException(LocalizedFormats.SCALE, scale);
128        }
129        this.scale = scale;
130        this.median = median;
131        solverAbsoluteAccuracy = inverseCumAccuracy;
132    }
133
134    /** {@inheritDoc} */
135    public double cumulativeProbability(double x) {
136        return 0.5 + (FastMath.atan((x - median) / scale) / FastMath.PI);
137    }
138
139    /**
140     * Access the median.
141     *
142     * @return the median for this distribution.
143     */
144    public double getMedian() {
145        return median;
146    }
147
148    /**
149     * Access the scale parameter.
150     *
151     * @return the scale parameter for this distribution.
152     */
153    public double getScale() {
154        return scale;
155    }
156
157    /** {@inheritDoc} */
158    public double density(double x) {
159        final double dev = x - median;
160        return (1 / FastMath.PI) * (scale / (dev * dev + scale * scale));
161    }
162
163    /**
164     * {@inheritDoc}
165     *
166     * Returns {@code Double.NEGATIVE_INFINITY} when {@code p == 0}
167     * and {@code Double.POSITIVE_INFINITY} when {@code p == 1}.
168     */
169    @Override
170    public double inverseCumulativeProbability(double p) throws OutOfRangeException {
171        double ret;
172        if (p < 0 || p > 1) {
173            throw new OutOfRangeException(p, 0, 1);
174        } else if (p == 0) {
175            ret = Double.NEGATIVE_INFINITY;
176        } else  if (p == 1) {
177            ret = Double.POSITIVE_INFINITY;
178        } else {
179            ret = median + scale * FastMath.tan(FastMath.PI * (p - .5));
180        }
181        return ret;
182    }
183
184    /** {@inheritDoc} */
185    @Override
186    protected double getSolverAbsoluteAccuracy() {
187        return solverAbsoluteAccuracy;
188    }
189
190    /**
191     * {@inheritDoc}
192     *
193     * The mean is always undefined no matter the parameters.
194     *
195     * @return mean (always Double.NaN)
196     */
197    public double getNumericalMean() {
198        return Double.NaN;
199    }
200
201    /**
202     * {@inheritDoc}
203     *
204     * The variance is always undefined no matter the parameters.
205     *
206     * @return variance (always Double.NaN)
207     */
208    public double getNumericalVariance() {
209        return Double.NaN;
210    }
211
212    /**
213     * {@inheritDoc}
214     *
215     * The lower bound of the support is always negative infinity no matter
216     * the parameters.
217     *
218     * @return lower bound of the support (always Double.NEGATIVE_INFINITY)
219     */
220    public double getSupportLowerBound() {
221        return Double.NEGATIVE_INFINITY;
222    }
223
224    /**
225     * {@inheritDoc}
226     *
227     * The upper bound of the support is always positive infinity no matter
228     * the parameters.
229     *
230     * @return upper bound of the support (always Double.POSITIVE_INFINITY)
231     */
232    public double getSupportUpperBound() {
233        return Double.POSITIVE_INFINITY;
234    }
235
236    /** {@inheritDoc} */
237    public boolean isSupportLowerBoundInclusive() {
238        return false;
239    }
240
241    /** {@inheritDoc} */
242    public boolean isSupportUpperBoundInclusive() {
243        return false;
244    }
245
246    /**
247     * {@inheritDoc}
248     *
249     * The support of this distribution is connected.
250     *
251     * @return {@code true}
252     */
253    public boolean isSupportConnected() {
254        return true;
255    }
256}