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 package org.hipparchus.distribution.discrete; 23 24 import java.io.Serializable; 25 26 import org.hipparchus.distribution.IntegerDistribution; 27 import org.hipparchus.exception.LocalizedCoreFormats; 28 import org.hipparchus.exception.MathIllegalArgumentException; 29 import org.hipparchus.exception.MathRuntimeException; 30 import org.hipparchus.util.FastMath; 31 import org.hipparchus.util.MathUtils; 32 33 /** 34 * Base class for integer-valued discrete distributions. 35 * <p> 36 * Default implementations are provided for some of the methods that 37 * do not vary from distribution to distribution. 38 */ 39 public abstract class AbstractIntegerDistribution implements IntegerDistribution, Serializable { 40 41 /** Serializable version identifier */ 42 private static final long serialVersionUID = 20160320L; 43 44 /** Empty constructor. 45 * <p> 46 * This constructor is not strictly necessary, but it prevents spurious 47 * javadoc warnings with JDK 18 and later. 48 * </p> 49 * @since 3.0 50 */ 51 public AbstractIntegerDistribution() { // NOPMD - unnecessary constructor added intentionally to make javadoc happy 52 // nothing to do 53 } 54 55 /** 56 * {@inheritDoc} 57 * 58 * The default implementation uses the identity 59 * <p> 60 * {@code P(x0 < X <= x1) = P(X <= x1) - P(X <= x0)} 61 */ 62 @Override 63 public double probability(int x0, int x1) throws MathIllegalArgumentException { 64 if (x1 < x0) { 65 throw new MathIllegalArgumentException(LocalizedCoreFormats.LOWER_ENDPOINT_ABOVE_UPPER_ENDPOINT, 66 x0, x1, true); 67 } 68 return cumulativeProbability(x1) - cumulativeProbability(x0); 69 } 70 71 /** 72 * {@inheritDoc} 73 * 74 * The default implementation returns 75 * <ul> 76 * <li>{@link #getSupportLowerBound()} for {@code p = 0},</li> 77 * <li>{@link #getSupportUpperBound()} for {@code p = 1}, and</li> 78 * <li>{@link #solveInverseCumulativeProbability(double, int, int)} for 79 * {@code 0 < p < 1}.</li> 80 * </ul> 81 */ 82 @Override 83 public int inverseCumulativeProbability(final double p) throws MathIllegalArgumentException { 84 MathUtils.checkRangeInclusive(p, 0, 1); 85 86 int lower = getSupportLowerBound(); 87 if (p == 0.0) { 88 return lower; 89 } 90 if (lower == Integer.MIN_VALUE) { 91 if (checkedCumulativeProbability(lower) >= p) { 92 return lower; 93 } 94 } else { 95 lower -= 1; // this ensures cumulativeProbability(lower) < p, which 96 // is important for the solving step 97 } 98 99 int upper = getSupportUpperBound(); 100 if (p == 1.0) { 101 return upper; 102 } 103 104 // use the one-sided Chebyshev inequality to narrow the bracket 105 // cf. AbstractRealDistribution.inverseCumulativeProbability(double) 106 final double mu = getNumericalMean(); 107 final double sigma = FastMath.sqrt(getNumericalVariance()); 108 final boolean chebyshevApplies = 109 !(Double.isInfinite(mu) || Double.isNaN(mu) || 110 Double.isInfinite(sigma) || Double.isNaN(sigma) || 111 sigma == 0.0); 112 if (chebyshevApplies) { 113 double k = FastMath.sqrt((1.0 - p) / p); 114 double tmp = mu - k * sigma; 115 if (tmp > lower) { 116 lower = ((int) FastMath.ceil(tmp)) - 1; 117 } 118 k = 1.0 / k; 119 tmp = mu + k * sigma; 120 if (tmp < upper) { 121 upper = ((int) FastMath.ceil(tmp)) - 1; 122 } 123 } 124 125 return solveInverseCumulativeProbability(p, lower, upper); 126 } 127 128 /** 129 * This is a utility function used by {@link 130 * #inverseCumulativeProbability(double)}. It assumes {@code 0 < p < 1} and 131 * that the inverse cumulative probability lies in the bracket {@code 132 * (lower, upper]}. The implementation does simple bisection to find the 133 * smallest {@code p}-quantile {@code inf{x in Z | P(X<=x) >= p}}. 134 * 135 * @param p the cumulative probability 136 * @param lower a value satisfying {@code cumulativeProbability(lower) < p} 137 * @param upper a value satisfying {@code p <= cumulativeProbability(upper)} 138 * @return the smallest {@code p}-quantile of this distribution 139 */ 140 protected int solveInverseCumulativeProbability(final double p, int lower, int upper) { 141 while (lower + 1 < upper) { 142 int xm = (lower + upper) / 2; 143 if (xm < lower || xm > upper) { 144 /* 145 * Overflow. 146 * There will never be an overflow in both calculation methods 147 * for xm at the same time 148 */ 149 xm = lower + (upper - lower) / 2; 150 } 151 152 double pm = checkedCumulativeProbability(xm); 153 if (pm >= p) { 154 upper = xm; 155 } else { 156 lower = xm; 157 } 158 } 159 return upper; 160 } 161 162 /** 163 * Computes the cumulative probability function and checks for {@code NaN} 164 * values returned. 165 * <p> 166 * Throws {@code MathRuntimeException} if the value is {@code NaN}. 167 * Rethrows any exception encountered evaluating the cumulative 168 * probability function. 169 * Throws {@code MathRuntimeException} if the cumulative 170 * probability function returns {@code NaN}. 171 * 172 * @param argument input value 173 * @return the cumulative probability 174 * @throws MathRuntimeException if the cumulative probability is {@code NaN} 175 */ 176 private double checkedCumulativeProbability(int argument) 177 throws MathRuntimeException { 178 double result = cumulativeProbability(argument); 179 if (Double.isNaN(result)) { 180 throw new MathRuntimeException(LocalizedCoreFormats.DISCRETE_CUMULATIVE_PROBABILITY_RETURNED_NAN, 181 argument); 182 } 183 return result; 184 } 185 186 /** 187 * {@inheritDoc} 188 * <p> 189 * The default implementation simply computes the logarithm of {@code probability(x)}. 190 */ 191 @Override 192 public double logProbability(int x) { 193 return FastMath.log(probability(x)); 194 } 195 }