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 package org.apache.commons.statistics.distribution;
18
19 import org.apache.commons.rng.UniformRandomProvider;
20 import org.apache.commons.rng.sampling.distribution.StableSampler;
21
22 /**
23 * Implementation of the Cauchy distribution.
24 *
25 * <p>The probability density function of \( X \) is:
26 *
27 * <p>\[ f(x; x_0, \gamma) = { 1 \over \pi \gamma } \left[ { \gamma^2 \over (x - x_0)^2 + \gamma^2 } \right] \]
28 *
29 * <p>for \( x_0 \) the location,
30 * \( \gamma > 0 \) the scale, and
31 * \( x \in (-\infty, \infty) \).
32 *
33 * @see <a href="https://en.wikipedia.org/wiki/Cauchy_distribution">Cauchy distribution (Wikipedia)</a>
34 * @see <a href="https://mathworld.wolfram.com/CauchyDistribution.html">Cauchy distribution (MathWorld)</a>
35 */
36 public final class CauchyDistribution extends AbstractContinuousDistribution {
37 /** The location of this distribution. */
38 private final double location;
39 /** The scale of this distribution. */
40 private final double scale;
41 /** Density factor (scale / pi). */
42 private final double scaleOverPi;
43 /** Density factor (scale^2). */
44 private final double scale2;
45
46 /**
47 * @param location Location parameter.
48 * @param scale Scale parameter.
49 */
50 private CauchyDistribution(double location,
51 double scale) {
52 this.scale = scale;
53 this.location = location;
54 scaleOverPi = scale / Math.PI;
55 scale2 = scale * scale;
56 }
57
58 /**
59 * Creates a Cauchy distribution.
60 *
61 * @param location Location parameter.
62 * @param scale Scale parameter.
63 * @return the distribution
64 * @throws IllegalArgumentException if {@code scale <= 0}.
65 */
66 public static CauchyDistribution of(double location,
67 double scale) {
68 if (scale <= 0) {
69 throw new DistributionException(DistributionException.NOT_STRICTLY_POSITIVE, scale);
70 }
71 return new CauchyDistribution(location, scale);
72 }
73
74 /**
75 * Gets the location parameter of this distribution.
76 *
77 * @return the location parameter.
78 */
79 public double getLocation() {
80 return location;
81 }
82
83 /**
84 * Gets the scale parameter of this distribution.
85 *
86 * @return the scale parameter.
87 */
88 public double getScale() {
89 return scale;
90 }
91
92 /** {@inheritDoc} */
93 @Override
94 public double density(double x) {
95 final double dev = x - location;
96 return scaleOverPi / (dev * dev + scale2);
97 }
98
99 /** {@inheritDoc} */
100 @Override
101 public double cumulativeProbability(double x) {
102 return cdf((x - location) / scale);
103 }
104
105 /** {@inheritDoc} */
106 @Override
107 public double survivalProbability(double x) {
108 return cdf(-(x - location) / scale);
109 }
110
111 /**
112 * Compute the CDF of the Cauchy distribution with location 0 and scale 1.
113 * @param x Point at which the CDF is evaluated
114 * @return CDF(x)
115 */
116 private static double cdf(double x) {
117 return 0.5 + (Math.atan(x) / Math.PI);
118 }
119
120 /**
121 * {@inheritDoc}
122 *
123 * <p>Returns {@link Double#NEGATIVE_INFINITY} when {@code p == 0}
124 * and {@link Double#POSITIVE_INFINITY} when {@code p == 1}.
125 */
126 @Override
127 public double inverseCumulativeProbability(double p) {
128 ArgumentUtils.checkProbability(p);
129 if (p == 0) {
130 return Double.NEGATIVE_INFINITY;
131 } else if (p == 1) {
132 return Double.POSITIVE_INFINITY;
133 }
134 return location + scale * Math.tan(Math.PI * (p - 0.5));
135 }
136
137 /**
138 * {@inheritDoc}
139 *
140 * <p>Returns {@link Double#NEGATIVE_INFINITY} when {@code p == 1}
141 * and {@link Double#POSITIVE_INFINITY} when {@code p == 0}.
142 */
143 @Override
144 public double inverseSurvivalProbability(double p) {
145 ArgumentUtils.checkProbability(p);
146 if (p == 1) {
147 return Double.NEGATIVE_INFINITY;
148 } else if (p == 0) {
149 return Double.POSITIVE_INFINITY;
150 }
151 return location - scale * Math.tan(Math.PI * (p - 0.5));
152 }
153
154 /**
155 * {@inheritDoc}
156 *
157 * <p>The mean is always undefined.
158 *
159 * @return {@link Double#NaN NaN}.
160 */
161 @Override
162 public double getMean() {
163 return Double.NaN;
164 }
165
166 /**
167 * {@inheritDoc}
168 *
169 * <p>The variance is always undefined.
170 *
171 * @return {@link Double#NaN NaN}.
172 */
173 @Override
174 public double getVariance() {
175 return Double.NaN;
176 }
177
178 /**
179 * {@inheritDoc}
180 *
181 * <p>The lower bound of the support is always negative infinity.
182 *
183 * @return {@linkplain Double#NEGATIVE_INFINITY negative infinity}.
184 */
185 @Override
186 public double getSupportLowerBound() {
187 return Double.NEGATIVE_INFINITY;
188 }
189
190 /**
191 * {@inheritDoc}
192 *
193 * <p>The upper bound of the support is always positive infinity.
194 *
195 * @return {@linkplain Double#POSITIVE_INFINITY positive infinity}.
196 */
197 @Override
198 public double getSupportUpperBound() {
199 return Double.POSITIVE_INFINITY;
200 }
201
202 /** {@inheritDoc} */
203 @Override
204 double getMedian() {
205 // Overridden for the probability(double, double) method.
206 // This is intentionally not a public method.
207 return location;
208 }
209
210 /** {@inheritDoc} */
211 @Override
212 public ContinuousDistribution.Sampler createSampler(final UniformRandomProvider rng) {
213 // Cauchy distribution =
214 // Stable distribution with alpha=1, beta=0, gamma=scale, delta=location
215 return StableSampler.of(rng, 1, 0, getScale(), getLocation())::sample;
216 }
217 }