1dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond/* 2dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Licensed to the Apache Software Foundation (ASF) under one or more 3dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * contributor license agreements. See the NOTICE file distributed with 4dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * this work for additional information regarding copyright ownership. 5dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * The ASF licenses this file to You under the Apache License, Version 2.0 6dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * (the "License"); you may not use this file except in compliance with 7dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * the License. You may obtain a copy of the License at 8dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 9dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * http://www.apache.org/licenses/LICENSE-2.0 10dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 11dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Unless required by applicable law or agreed to in writing, software 12dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * distributed under the License is distributed on an "AS IS" BASIS, 13dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 14dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * See the License for the specific language governing permissions and 15dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * limitations under the License. 16dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 17dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 18dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondpackage org.apache.commons.math.distribution; 19dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 20dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport java.io.Serializable; 21dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 22dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.MathRuntimeException; 23dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.exception.util.LocalizedFormats; 24dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.special.Gamma; 25dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.util.FastMath; 26dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 27dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond/** 28dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Default implementation of 29dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * {@link org.apache.commons.math.distribution.WeibullDistribution}. 30dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 31dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @since 1.1 32dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @version $Revision: 1054524 $ $Date: 2011-01-03 05:59:18 +0100 (lun. 03 janv. 2011) $ 33dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 34dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondpublic class WeibullDistributionImpl extends AbstractContinuousDistribution 35dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond implements WeibullDistribution, Serializable { 36dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 37dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 38dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Default inverse cumulative probability accuracy 39dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @since 2.1 40dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 41dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9; 42dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 43dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Serializable version identifier */ 44dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private static final long serialVersionUID = 8589540077390120676L; 45dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 46dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** The shape parameter. */ 47dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private double shape; 48dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 49dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** The scale parameter. */ 50dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private double scale; 51dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 52dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Inverse cumulative probability accuracy */ 53dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private final double solverAbsoluteAccuracy; 54dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 55dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Cached numerical mean */ 56dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private double numericalMean = Double.NaN; 57dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 58dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Whether or not the numerical mean has been calculated */ 59dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private boolean numericalMeanIsCalculated = false; 60dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 61dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Cached numerical variance */ 62dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private double numericalVariance = Double.NaN; 63dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 64dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Whether or not the numerical variance has been calculated */ 65dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private boolean numericalVarianceIsCalculated = false; 66dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 67dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 68dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Creates weibull distribution with the given shape and scale and a 69dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * location equal to zero. 70dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param alpha the shape parameter. 71dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param beta the scale parameter. 72dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 73dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public WeibullDistributionImpl(double alpha, double beta){ 74dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this(alpha, beta, DEFAULT_INVERSE_ABSOLUTE_ACCURACY); 75dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 76dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 77dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 78dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Creates weibull distribution with the given shape, scale and inverse 79dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * cumulative probability accuracy and a location equal to zero. 80dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param alpha the shape parameter. 81dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param beta the scale parameter. 82dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param inverseCumAccuracy the maximum absolute error in inverse cumulative probability estimates 83dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * (defaults to {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}) 84dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @since 2.1 85dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 86dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public WeibullDistributionImpl(double alpha, double beta, double inverseCumAccuracy){ 87dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond super(); 88dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond setShapeInternal(alpha); 89dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond setScaleInternal(beta); 90dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond solverAbsoluteAccuracy = inverseCumAccuracy; 91dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 92dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 93dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 94dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * For this distribution, X, this method returns P(X < <code>x</code>). 95dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param x the value at which the CDF is evaluated. 96dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return CDF evaluated at <code>x</code>. 97dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 98dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public double cumulativeProbability(double x) { 99dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double ret; 100dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (x <= 0.0) { 101dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond ret = 0.0; 102dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } else { 103dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond ret = 1.0 - FastMath.exp(-FastMath.pow(x / scale, shape)); 104dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 105dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return ret; 106dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 107dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 108dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 109dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Access the shape parameter. 110dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return the shape parameter. 111dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 112dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public double getShape() { 113dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return shape; 114dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 115dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 116dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 117dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Access the scale parameter. 118dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return the scale parameter. 119dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 120dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public double getScale() { 121dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return scale; 122dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 123dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 124dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 125dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Returns the probability density for a particular point. 126dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 127dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param x The point at which the density should be computed. 128dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return The pdf at point x. 129dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @since 2.1 130dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 131dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond @Override 132dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public double density(double x) { 133dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (x < 0) { 134dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return 0; 135dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 136dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 137dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double xscale = x / scale; 138dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double xscalepow = FastMath.pow(xscale, shape - 1); 139dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 140dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /* 141dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * FastMath.pow(x / scale, shape) = 142dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * FastMath.pow(xscale, shape) = 143dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * FastMath.pow(xscale, shape - 1) * xscale 144dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 145dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double xscalepowshape = xscalepow * xscale; 146dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 147dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return (shape / scale) * xscalepow * FastMath.exp(-xscalepowshape); 148dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 149dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 150dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 151dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * For this distribution, X, this method returns the critical point x, such 152dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * that P(X < x) = <code>p</code>. 153dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p> 154dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Returns <code>Double.NEGATIVE_INFINITY</code> for p=0 and 155dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <code>Double.POSITIVE_INFINITY</code> for p=1.</p> 156dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 157dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param p the desired probability 158dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return x, such that P(X < x) = <code>p</code> 159dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @throws IllegalArgumentException if <code>p</code> is not a valid 160dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * probability. 161dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 162dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond @Override 163dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public double inverseCumulativeProbability(double p) { 164dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double ret; 165dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (p < 0.0 || p > 1.0) { 166dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw MathRuntimeException.createIllegalArgumentException( 167dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond LocalizedFormats.OUT_OF_RANGE_SIMPLE, p, 0.0, 1.0); 168dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } else if (p == 0) { 169dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond ret = 0.0; 170dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } else if (p == 1) { 171dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond ret = Double.POSITIVE_INFINITY; 172dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } else { 173dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond ret = scale * FastMath.pow(-FastMath.log(1.0 - p), 1.0 / shape); 174dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 175dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return ret; 176dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 177dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 178dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 179dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Modify the shape parameter. 180dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param alpha the new shape parameter value. 181dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @deprecated as of 2.1 (class will become immutable in 3.0) 182dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 183dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond @Deprecated 184dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public void setShape(double alpha) { 185dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond setShapeInternal(alpha); 186dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond invalidateParameterDependentMoments(); 187dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 188dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 189dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Modify the shape parameter. 190dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param alpha the new shape parameter value. 191dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 192dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private void setShapeInternal(double alpha) { 193dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (alpha <= 0.0) { 194dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw MathRuntimeException.createIllegalArgumentException( 195dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond LocalizedFormats.NOT_POSITIVE_SHAPE, 196dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond alpha); 197dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 198dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.shape = alpha; 199dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 200dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 201dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 202dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Modify the scale parameter. 203dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param beta the new scale parameter value. 204dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @deprecated as of 2.1 (class will become immutable in 3.0) 205dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 206dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond @Deprecated 207dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public void setScale(double beta) { 208dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond setScaleInternal(beta); 209dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond invalidateParameterDependentMoments(); 210dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 211dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 212dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Modify the scale parameter. 213dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param beta the new scale parameter value. 214dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 215dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private void setScaleInternal(double beta) { 216dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (beta <= 0.0) { 217dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw MathRuntimeException.createIllegalArgumentException( 218dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond LocalizedFormats.NOT_POSITIVE_SCALE, 219dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond beta); 220dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 221dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.scale = beta; 222dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 223dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 224dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 225dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Access the domain value lower bound, based on <code>p</code>, used to 226dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * bracket a CDF root. This method is used by 227dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * {@link #inverseCumulativeProbability(double)} to find critical values. 228dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 229dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param p the desired probability for the critical value 230dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return domain value lower bound, i.e. 231dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * P(X < <i>lower bound</i>) < <code>p</code> 232dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 233dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond @Override 234dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected double getDomainLowerBound(double p) { 235dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return 0.0; 236dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 237dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 238dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 239dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Access the domain value upper bound, based on <code>p</code>, used to 240dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * bracket a CDF root. This method is used by 241dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * {@link #inverseCumulativeProbability(double)} to find critical values. 242dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 243dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param p the desired probability for the critical value 244dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return domain value upper bound, i.e. 245dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * P(X < <i>upper bound</i>) > <code>p</code> 246dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 247dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond @Override 248dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected double getDomainUpperBound(double p) { 249dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return Double.MAX_VALUE; 250dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 251dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 252dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 253dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Access the initial domain value, based on <code>p</code>, used to 254dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * bracket a CDF root. This method is used by 255dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * {@link #inverseCumulativeProbability(double)} to find critical values. 256dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 257dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param p the desired probability for the critical value 258dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return initial domain value 259dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 260dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond @Override 261dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected double getInitialDomain(double p) { 262dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // use median 263dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return FastMath.pow(scale * FastMath.log(2.0), 1.0 / shape); 264dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 265dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 266dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 267dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Return the absolute accuracy setting of the solver used to estimate 268dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * inverse cumulative probabilities. 269dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 270dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return the solver absolute accuracy 271dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @since 2.1 272dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 273dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond @Override 274dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected double getSolverAbsoluteAccuracy() { 275dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return solverAbsoluteAccuracy; 276dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 277dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 278dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 279dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Returns the lower bound of the support for the distribution. 280dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 281dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * The lower bound of the support is always 0 no matter the parameters. 282dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 283dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return lower bound of the support (always 0) 284dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @since 2.2 285dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 286dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public double getSupportLowerBound() { 287dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return 0; 288dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 289dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 290dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 291dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Returns the upper bound of the support for the distribution. 292dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 293dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * The upper bound of the support is always positive infinity 294dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * no matter the parameters. 295dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 296dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return upper bound of the support (always Double.POSITIVE_INFINITY) 297dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @since 2.2 298dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 299dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public double getSupportUpperBound() { 300dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return Double.POSITIVE_INFINITY; 301dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 302dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 303dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 304dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Calculates the mean. 305dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 306dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * The mean is <code>scale * Gamma(1 + (1 / shape))</code> 307dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * where <code>Gamma(...)</code> is the Gamma-function 308dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 309dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return the mean 310dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @since 2.2 311dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 312dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond protected double calculateNumericalMean() { 313dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double sh = getShape(); 314dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double sc = getScale(); 315dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 316dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return sc * FastMath.exp(Gamma.logGamma(1 + (1 / sh))); 317dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 318dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 319dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 320dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Calculates the variance. 321dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 322dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * The variance is 323dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <code>scale^2 * Gamma(1 + (2 / shape)) - mean^2</code> 324dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * where <code>Gamma(...)</code> is the Gamma-function 325dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 326dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return the variance 327dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @since 2.2 328dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 329dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private double calculateNumericalVariance() { 330dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double sh = getShape(); 331dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double sc = getScale(); 332dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double mn = getNumericalMean(); 333dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 334dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return (sc * sc) * 335dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond FastMath.exp(Gamma.logGamma(1 + (2 / sh))) - 336dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond (mn * mn); 337dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 338dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 339dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 340dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Returns the mean of the distribution. 341dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 342dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return the mean or Double.NaN if it's not defined 343dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @since 2.2 344dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 345dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public double getNumericalMean() { 346dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (!numericalMeanIsCalculated) { 347dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond numericalMean = calculateNumericalMean(); 348dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond numericalMeanIsCalculated = true; 349dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 350dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 351dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return numericalMean; 352dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 353dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 354dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 355dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Returns the variance of the distribution. 356dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 357dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return the variance (possibly Double.POSITIVE_INFINITY as 358dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * for certain cases in {@link TDistributionImpl}) or 359dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Double.NaN if it's not defined 360dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @since 2.2 361dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 362dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public double getNumericalVariance() { 363dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (!numericalVarianceIsCalculated) { 364dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond numericalVariance = calculateNumericalVariance(); 365dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond numericalVarianceIsCalculated = true; 366dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 367dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 368dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return numericalVariance; 369dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 370dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 371dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 372dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Invalidates the cached mean and variance. 373dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 374dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private void invalidateParameterDependentMoments() { 375dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond numericalMeanIsCalculated = false; 376dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond numericalVarianceIsCalculated = false; 377dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 378dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond} 379