View Javadoc
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 }