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 * https://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
18 /*
19 * This is not the original file distributed by the Apache Software Foundation
20 * It has been modified by the Hipparchus project
21 */
22
23 package org.hipparchus.distribution.continuous;
24
25 import org.hipparchus.exception.LocalizedCoreFormats;
26 import org.hipparchus.exception.MathIllegalArgumentException;
27 import org.hipparchus.util.FastMath;
28
29 /**
30 * Implementation of the Pareto distribution.
31 * <p>
32 * <strong>Parameters:</strong>
33 * The probability distribution function of {@code X} is given by (for {@code x >= k}):
34 * </p>
35 * <pre>
36 * α * k^α / x^(α + 1)
37 * </pre>
38 * <ul>
39 * <li>{@code k} is the <em>scale</em> parameter: this is the minimum possible value of {@code X},</li>
40 * <li>{@code α} is the <em>shape</em> parameter: this is the Pareto index</li>
41 * </ul>
42 *
43 * @see <a href="http://en.wikipedia.org/wiki/Pareto_distribution">
44 * Pareto distribution (Wikipedia)</a>
45 * @see <a href="http://mathworld.wolfram.com/ParetoDistribution.html">
46 * Pareto distribution (MathWorld)</a>
47 */
48 public class ParetoDistribution extends AbstractRealDistribution {
49
50 /** Serializable version identifier. */
51 private static final long serialVersionUID = 20130424L;
52
53 /** The scale parameter of this distribution. */
54 private final double scale;
55 /** The shape parameter of this distribution. */
56 private final double shape;
57
58 /**
59 * Create a Pareto distribution with a scale of {@code 1} and a shape of {@code 1}.
60 */
61 public ParetoDistribution() {
62 this(1, 1);
63 }
64
65 /**
66 * Create a Pareto distribution using the specified scale and shape.
67 *
68 * @param scale the scale parameter of this distribution
69 * @param shape the shape parameter of this distribution
70 * @throws MathIllegalArgumentException if {@code scale <= 0} or {@code shape <= 0}.
71 */
72 public ParetoDistribution(double scale, double shape)
73 throws MathIllegalArgumentException {
74 this(scale, shape, DEFAULT_SOLVER_ABSOLUTE_ACCURACY);
75 }
76
77 /**
78 * Creates a Pareto distribution.
79 *
80 * @param scale Scale parameter of this distribution.
81 * @param shape Shape parameter of this distribution.
82 * @param inverseCumAccuracy Inverse cumulative probability accuracy.
83 * @throws MathIllegalArgumentException if {@code scale <= 0} or {@code shape <= 0}.
84 */
85 public ParetoDistribution(double scale,
86 double shape,
87 double inverseCumAccuracy)
88 throws MathIllegalArgumentException {
89 super(inverseCumAccuracy);
90
91 if (scale <= 0) {
92 throw new MathIllegalArgumentException(LocalizedCoreFormats.SCALE, scale);
93 }
94
95 if (shape <= 0) {
96 throw new MathIllegalArgumentException(LocalizedCoreFormats.SHAPE, shape);
97 }
98
99 this.scale = scale;
100 this.shape = shape;
101 }
102
103 /**
104 * Returns the scale parameter of this distribution.
105 *
106 * @return the scale parameter
107 */
108 public double getScale() {
109 return scale;
110 }
111
112 /**
113 * Returns the shape parameter of this distribution.
114 *
115 * @return the shape parameter
116 */
117 public double getShape() {
118 return shape;
119 }
120
121 /**
122 * {@inheritDoc}
123 * <p>
124 * For scale {@code k}, and shape {@code α} of this distribution, the PDF
125 * is given by
126 * <ul>
127 * <li>{@code 0} if {@code x < k},</li>
128 * <li>{@code α * k^α / x^(α + 1)} otherwise.</li>
129 * </ul>
130 */
131 @Override
132 public double density(double x) {
133 if (x < scale) {
134 return 0;
135 }
136 return FastMath.pow(scale, shape) / FastMath.pow(x, shape + 1) * shape;
137 }
138
139 /** {@inheritDoc}
140 *
141 * See documentation of {@link #density(double)} for computation details.
142 */
143 @Override
144 public double logDensity(double x) {
145 if (x < scale) {
146 return Double.NEGATIVE_INFINITY;
147 }
148 return FastMath.log(scale) * shape - FastMath.log(x) * (shape + 1) + FastMath.log(shape);
149 }
150
151 /**
152 * {@inheritDoc}
153 * <p>
154 * For scale {@code k}, and shape {@code α} of this distribution, the CDF is given by
155 * <ul>
156 * <li>{@code 0} if {@code x < k},</li>
157 * <li>{@code 1 - (k / x)^α} otherwise.</li>
158 * </ul>
159 */
160 @Override
161 public double cumulativeProbability(double x) {
162 if (x <= scale) {
163 return 0;
164 }
165 return 1 - FastMath.pow(scale / x, shape);
166 }
167
168 /**
169 * {@inheritDoc}
170 * <p>
171 * For scale {@code k} and shape {@code α}, the mean is given by
172 * <ul>
173 * <li>{@code ∞} if {@code α <= 1},</li>
174 * <li>{@code α * k / (α - 1)} otherwise.</li>
175 * </ul>
176 */
177 @Override
178 public double getNumericalMean() {
179 if (shape <= 1) {
180 return Double.POSITIVE_INFINITY;
181 }
182 return shape * scale / (shape - 1);
183 }
184
185 /**
186 * {@inheritDoc}
187 * <p>
188 * For scale {@code k} and shape {@code α}, the variance is given by
189 * <ul>
190 * <li>{@code ∞} if {@code 1 < α <= 2},</li>
191 * <li>{@code k^2 * α / ((α - 1)^2 * (α - 2))} otherwise.</li>
192 * </ul>
193 */
194 @Override
195 public double getNumericalVariance() {
196 if (shape <= 2) {
197 return Double.POSITIVE_INFINITY;
198 }
199 double s = shape - 1;
200 return scale * scale * shape / (s * s) / (shape - 2);
201 }
202
203 /**
204 * {@inheritDoc}
205 * <p>
206 * The lower bound of the support is equal to the scale parameter {@code k}.
207 *
208 * @return lower bound of the support
209 */
210 @Override
211 public double getSupportLowerBound() {
212 return scale;
213 }
214
215 /**
216 * {@inheritDoc}
217 * <p>
218 * The upper bound of the support is always positive infinity no matter the parameters.
219 *
220 * @return upper bound of the support (always {@code Double.POSITIVE_INFINITY})
221 */
222 @Override
223 public double getSupportUpperBound() {
224 return Double.POSITIVE_INFINITY;
225 }
226
227 /**
228 * {@inheritDoc}
229 * <p>
230 * The support of this distribution is connected.
231 *
232 * @return {@code true}
233 */
234 @Override
235 public boolean isSupportConnected() {
236 return true;
237 }
238 }