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    package org.apache.commons.math.distribution;
018    
019    import java.io.Serializable;
020    
021    import org.apache.commons.math.MathException;
022    import org.apache.commons.math.MathRuntimeException;
023    
024    
025    /**
026     * Base class for integer-valued discrete distributions.  Default
027     * implementations are provided for some of the methods that do not vary
028     * from distribution to distribution.
029     *  
030     * @version $Revision: 772119 $ $Date: 2009-05-06 05:43:28 -0400 (Wed, 06 May 2009) $
031     */
032    public abstract class AbstractIntegerDistribution extends AbstractDistribution
033        implements IntegerDistribution, Serializable {
034            
035        /** Serializable version identifier */
036        private static final long serialVersionUID = -1146319659338487221L;
037        
038        /**
039         * Default constructor.
040         */
041        protected AbstractIntegerDistribution() {
042            super();
043        }
044        
045        /**
046         * For a random variable X whose values are distributed according
047         * to this distribution, this method returns P(X ≤ x).  In other words,
048         * this method represents the  (cumulative) distribution function, or
049         * CDF, for this distribution.
050         * <p>
051         * If <code>x</code> does not represent an integer value, the CDF is 
052         * evaluated at the greatest integer less than x.
053         * 
054         * @param x the value at which the distribution function is evaluated.
055         * @return cumulative probability that a random variable with this
056         * distribution takes a value less than or equal to <code>x</code>
057         * @throws MathException if the cumulative probability can not be
058         * computed due to convergence or other numerical errors.
059         */
060        public double cumulativeProbability(double x) throws MathException {
061            return cumulativeProbability((int) Math.floor(x));  
062        }
063        
064        /**
065         * For a random variable X whose values are distributed according
066         * to this distribution, this method returns P(x0 &le; X &le; x1).
067         * 
068         * @param x0 the (inclusive) lower bound
069         * @param x1 the (inclusive) upper bound
070         * @return the probability that a random variable with this distribution
071         * will take a value between <code>x0</code> and <code>x1</code>,
072         * including the endpoints.
073         * @throws MathException if the cumulative probability can not be
074         * computed due to convergence or other numerical errors.
075         * @throws IllegalArgumentException if <code>x0 > x1</code>
076         */
077        @Override
078        public double cumulativeProbability(double x0, double x1)
079            throws MathException {
080            if (x0 > x1) {
081                throw MathRuntimeException.createIllegalArgumentException(
082                      "lower endpoint ({0}) must be less than or equal to upper endpoint ({1})",
083                      x0, x1);
084            }
085            if (Math.floor(x0) < x0) {
086                return cumulativeProbability(((int) Math.floor(x0)) + 1,
087                   (int) Math.floor(x1)); // don't want to count mass below x0
088            } else { // x0 is mathematical integer, so use as is
089                return cumulativeProbability((int) Math.floor(x0),
090                    (int) Math.floor(x1)); 
091            }
092        }
093        
094        /**
095         * For a random variable X whose values are distributed according
096         * to this distribution, this method returns P(X &le; x).  In other words,
097         * this method represents the probability distribution function, or PDF,
098         * for this distribution.
099         * 
100         * @param x the value at which the PDF is evaluated.
101         * @return PDF for this distribution. 
102         * @throws MathException if the cumulative probability can not be
103         *            computed due to convergence or other numerical errors.
104         */
105        abstract public double cumulativeProbability(int x) throws MathException;
106        
107        /**
108         * For a random variable X whose values are distributed according
109         * to this distribution, this method returns P(X = x). In other words, this
110         * method represents the probability mass function,  or PMF, for the distribution.
111         * <p>
112         * If <code>x</code> does not represent an integer value, 0 is returned.
113         * 
114         * @param x the value at which the probability density function is evaluated
115         * @return the value of the probability density function at x
116         */
117        public double probability(double x) {
118            double fl = Math.floor(x);
119            if (fl == x) {
120                return this.probability((int) x);
121            } else {
122                return 0;
123            }
124        }
125        
126        /**
127        * For a random variable X whose values are distributed according
128         * to this distribution, this method returns P(x0 &le; X &le; x1).
129         * 
130         * @param x0 the inclusive, lower bound
131         * @param x1 the inclusive, upper bound
132         * @return the cumulative probability. 
133         * @throws MathException if the cumulative probability can not be
134         *            computed due to convergence or other numerical errors.
135         * @throws IllegalArgumentException if x0 > x1
136         */
137        public double cumulativeProbability(int x0, int x1) throws MathException {
138            if (x0 > x1) {
139                throw MathRuntimeException.createIllegalArgumentException(
140                      "lower endpoint ({0}) must be less than or equal to upper endpoint ({1})",
141                      x0, x1);
142            }
143            return cumulativeProbability(x1) - cumulativeProbability(x0 - 1);
144        }
145        
146        /**
147         * For a random variable X whose values are distributed according
148         * to this distribution, this method returns the largest x, such
149         * that P(X &le; x) &le; <code>p</code>.
150         *
151         * @param p the desired probability
152         * @return the largest x such that P(X &le; x) <= p
153         * @throws MathException if the inverse cumulative probability can not be
154         *            computed due to convergence or other numerical errors.
155         * @throws IllegalArgumentException if p < 0 or p > 1
156         */
157        public int inverseCumulativeProbability(final double p) throws MathException{
158            if (p < 0.0 || p > 1.0) {
159                throw MathRuntimeException.createIllegalArgumentException(
160                      "{0} out of [{1}, {2}] range", p, 0.0, 1.0);
161            }
162            
163            // by default, do simple bisection.
164            // subclasses can override if there is a better method.
165            int x0 = getDomainLowerBound(p);
166            int x1 = getDomainUpperBound(p);
167            double pm;
168            while (x0 < x1) {
169                int xm = x0 + (x1 - x0) / 2;
170                pm = cumulativeProbability(xm);
171                if (pm > p) {
172                    // update x1
173                    if (xm == x1) {
174                        // this can happen with integer division
175                        // simply decrement x1
176                        --x1;
177                    } else {
178                        // update x1 normally
179                        x1 = xm;
180                    }
181                } else {
182                    // update x0
183                    if (xm == x0) {
184                        // this can happen with integer division
185                        // simply increment x0
186                        ++x0;
187                    } else {
188                        // update x0 normally
189                        x0 = xm;
190                    }
191                }
192            }
193            
194            // insure x0 is the correct critical point
195            pm = cumulativeProbability(x0);
196            while (pm > p) {
197                --x0;
198                pm = cumulativeProbability(x0);
199            }
200        
201            return x0;        
202        }
203        
204        /**
205         * Access the domain value lower bound, based on <code>p</code>, used to
206         * bracket a PDF root.  This method is used by
207         * {@link #inverseCumulativeProbability(double)} to find critical values.
208         * 
209         * @param p the desired probability for the critical value
210         * @return domain value lower bound, i.e.
211         *         P(X &lt; <i>lower bound</i>) &lt; <code>p</code> 
212         */
213        protected abstract int getDomainLowerBound(double p);
214        
215        /**
216         * Access the domain value upper bound, based on <code>p</code>, used to
217         * bracket a PDF root.  This method is used by
218         * {@link #inverseCumulativeProbability(double)} to find critical values.
219         * 
220         * @param p the desired probability for the critical value
221         * @return domain value upper bound, i.e.
222         *         P(X &lt; <i>upper bound</i>) &gt; <code>p</code> 
223         */
224        protected abstract int getDomainUpperBound(double p);
225    }