1 // -*- C++ -*-
2 /**
3 * @brief The InverseChiSquare distribution
4 *
5 * Copyright 2005-2021 Airbus-EDF-IMACS-ONERA-Phimeca
6 *
7 * This library is free software: you can redistribute it and/or modify
8 * it under the terms of the GNU Lesser General Public License as published by
9 * the Free Software Foundation, either version 3 of the License, or
10 * (at your option) any later version.
11 *
12 * This library is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU Lesser General Public License for more details.
16 *
17 * You should have received a copy of the GNU Lesser General Public License
18 * along with this library. If not, see <http://www.gnu.org/licenses/>.
19 *
20 */
21 #include <cmath>
22 #include "openturns/InverseChiSquare.hxx"
23 #include "openturns/RandomGenerator.hxx"
24 #include "openturns/SpecFunc.hxx"
25 #include "openturns/DistFunc.hxx"
26 #include "openturns/PersistentObjectFactory.hxx"
27 #include "openturns/Distribution.hxx"
28
29 BEGIN_NAMESPACE_OPENTURNS
30
31 CLASSNAMEINIT(InverseChiSquare)
32
33 static const Factory<InverseChiSquare> Factory_InverseChiSquare;
34
35 /* Default constructor */
InverseChiSquare()36 InverseChiSquare::InverseChiSquare()
37 : ContinuousDistribution()
38 , nu_(1.0)
39 , normalizationFactor_(0.0)
40 {
41 setName("InverseChiSquare");
42 setDimension(1);
43 computeRange();
44 }
45
46 /* Parameters constructor */
InverseChiSquare(const Scalar nu)47 InverseChiSquare::InverseChiSquare(const Scalar nu)
48 : ContinuousDistribution()
49 , nu_(0.0)
50 , normalizationFactor_(0.0)
51 {
52 setName("InverseChiSquare");
53 setNu(nu);
54 setDimension(1);
55 }
56
57 /* Comparison operator */
operator ==(const InverseChiSquare & other) const58 Bool InverseChiSquare::operator ==(const InverseChiSquare & other) const
59 {
60 if (this == &other) return true;
61 return nu_ == other.nu_;
62 }
63
equals(const DistributionImplementation & other) const64 Bool InverseChiSquare::equals(const DistributionImplementation & other) const
65 {
66 const InverseChiSquare* p_other = dynamic_cast<const InverseChiSquare*>(&other);
67 return p_other && (*this == *p_other);
68 }
69
70 /* String converter */
__repr__() const71 String InverseChiSquare::__repr__() const
72 {
73 OSS oss;
74 oss << "class=" << InverseChiSquare::GetClassName()
75 << " name=" << getName()
76 << " dimension=" << getDimension()
77 << " nu=" << nu_;
78 return oss;
79 }
80
__str__(const String &) const81 String InverseChiSquare::__str__(const String & ) const
82 {
83 OSS oss;
84 oss << getClassName() << "(nu = " << nu_ << ")";
85 return oss;
86 }
87
88 /* K accessor */
setNu(const Scalar nu)89 void InverseChiSquare::setNu(const Scalar nu)
90 {
91 if (!(nu > 0.0)) throw InvalidArgumentException(HERE) << "Nu MUST be positive";
92 if (nu != nu_)
93 {
94 nu_ = nu;
95 computeRange();
96 update();
97 }
98 }
99
getNu() const100 Scalar InverseChiSquare::getNu() const
101 {
102 return nu_;
103 }
104
105 /* Virtual constructor */
clone() const106 InverseChiSquare * InverseChiSquare::clone() const
107 {
108 return new InverseChiSquare(*this);
109 }
110
111 /* Compute the numerical range of the distribution given the parameters values */
computeRange()112 void InverseChiSquare::computeRange()
113 {
114 const Point lowerBound(1, 0.0);
115 const Point upperBound(computeUpperBound());
116 const Interval::BoolCollection finiteLowerBound(1, true);
117 const Interval::BoolCollection finiteUpperBound(1, false);
118 setRange(Interval(lowerBound, upperBound, finiteLowerBound, finiteUpperBound));
119 }
120
121 /** Update the derivative attributes */
update()122 void InverseChiSquare::update()
123 {
124 // For large k we use the following normalization factor:
125 // normalizationFactor = log(2*k^{k-1}/Gamma(k))
126 // = log(2) + (k+1)log(k) - log(Gamma(k))
127 // which is expanded wrt k
128 const Scalar k = 0.5 * nu_;
129 if (k >= 6.9707081224932495879)
130 {
131 static const Scalar alpha[10] = {0.91893853320467274177, 0.83333333333333333333e-1, -0.27777777777777777778e-2, 0.79365079365079365079e-3, -0.59523809523809523810e-3, 0.84175084175084175084e-3, -0.19175269175269175269e-2, 0.64102564102564102564e-2, -0.29550653594771241830e-1, 0.17964437236883057316};
132 const Scalar ik = 1.0 / k;
133 const Scalar ik2 = ik * ik;
134 normalizationFactor_ = std::log(2.0) + k + 1.5 * std::log(k) - (alpha[0] + ik * (alpha[1] + ik2 * (alpha[2] + ik2 * (alpha[3] + ik2 * (alpha[4] + ik2 * (alpha[5] + ik2 * (alpha[6] + ik2 * (alpha[7] + ik2 * (alpha[8] + ik2 * alpha[9])))))))));
135 }
136 // For small k, the normalization factor is:
137 // normalizationFactor = log(2/Gamma(k))
138 // = log(2) - log(Gamma(k))
139 else normalizationFactor_ = std::log(2.0) - SpecFunc::LnGamma(k);
140 isAlreadyComputedMean_ = false;
141 isAlreadyComputedCovariance_ = false;
142 }
143
144
145 /* Get one realization of the distribution */
getRealization() const146 Point InverseChiSquare::getRealization() const
147 {
148 return Point(1, 1.0 / (2.0 * DistFunc::rGamma(0.5 * nu_)));
149 }
150
151
152 /* Get the DDF of the distribution */
computeDDF(const Point & point) const153 Point InverseChiSquare::computeDDF(const Point & point) const
154 {
155 if (point.getDimension() != 1) throw InvalidArgumentException(HERE) << "Error: the given point must have dimension=1, here dimension=" << point.getDimension();
156
157 const Scalar x = point[0];
158 if (x <= 0.0) return Point(1, 0.0);
159 return Point(1, (1.0 / (2.0 * x) - (0.5 * nu_ + 1.0)) * computePDF(point) / x);
160 }
161
162
163 /* Get the PDF of the distribution */
computePDF(const Point & point) const164 Scalar InverseChiSquare::computePDF(const Point & point) const
165 {
166 if (point.getDimension() != 1) throw InvalidArgumentException(HERE) << "Error: the given point must have dimension=1, here dimension=" << point.getDimension();
167
168 if (point[0] <= 0.0) return 0.0;
169 return std::exp(computeLogPDF(point));
170 }
171
computeLogPDF(const Point & point) const172 Scalar InverseChiSquare::computeLogPDF(const Point & point) const
173 {
174 if (point.getDimension() != 1) throw InvalidArgumentException(HERE) << "Error: the given point must have dimension=1, here dimension=" << point.getDimension();
175
176 // From textbook, we have log(PDF(x)) = log(lambda)-log(Gamma(k))-(k+1)*log(lambda*x)-1/(lambda*x)
177 const Scalar u = 2.0 * point[0];
178 if (u <= 0.0) return SpecFunc::LowestScalar;
179 // Use asymptotic expansion for large k
180 // Here log(PDF(x)) = L - (k-1)*log(k)-(k+1)*log(2*x)-1/(2*x)
181 const Scalar k = 0.5 * nu_;
182 if (k >= 6.9707081224932495879) return normalizationFactor_ - (k + 1.0) * std::log(k * u) - 1.0 / u;
183 return normalizationFactor_ - (k + 1.0) * std::log(u) - 1.0 / u;
184 }
185
186 /* Get the CDF of the distribution */
computeCDF(const Point & point) const187 Scalar InverseChiSquare::computeCDF(const Point & point) const
188 {
189 if (point.getDimension() != 1) throw InvalidArgumentException(HERE) << "Error: the given point must have dimension=1, here dimension=" << point.getDimension();
190
191 const Scalar x = point[0];
192 // No test here as the CDF is continuous for all k
193 if (x <= 0.0) return 0.0;
194 return DistFunc::pGamma(0.5 * nu_, 0.5 / x, true);
195 }
196
computeComplementaryCDF(const Point & point) const197 Scalar InverseChiSquare::computeComplementaryCDF(const Point & point) const
198 {
199 if (point.getDimension() != 1) throw InvalidArgumentException(HERE) << "Error: the given point must have dimension=1, here dimension=" << point.getDimension();
200
201 const Scalar x = point[0];
202 // No test here as the CDF is continuous for all k
203 if (x <= 0.0) return 1.0;
204 return DistFunc::pGamma(0.5 * nu_, 0.5 / x);
205 }
206
207 /* Compute the entropy of the distribution */
computeEntropy() const208 Scalar InverseChiSquare::computeEntropy() const
209 {
210 return 0.5 * nu_ - M_LN2 + SpecFunc::LogGamma(0.5 * nu_) - (1.0 + 0.5 * nu_) * SpecFunc::Psi(0.5 * nu_);
211 }
212
213 /* Get the characteristic function of the distribution, i.e. phi(u) = E(exp(I*u*X)) */
computeCharacteristicFunction(const Scalar x) const214 Complex InverseChiSquare::computeCharacteristicFunction(const Scalar x) const
215 {
216 return DistributionImplementation::computeCharacteristicFunction(x);
217 }
218
computeLogCharacteristicFunction(const Scalar x) const219 Complex InverseChiSquare::computeLogCharacteristicFunction(const Scalar x) const
220 {
221 return DistributionImplementation::computeLogCharacteristicFunction(x);
222 }
223
224 /* Get the PDFGradient of the distribution */
computePDFGradient(const Point & point) const225 Point InverseChiSquare::computePDFGradient(const Point & point) const
226 {
227 if (point.getDimension() != 1) throw InvalidArgumentException(HERE) << "Error: the given point must have dimension=1, here dimension=" << point.getDimension();
228
229 Point pdfGradient(2);
230 const Scalar x = point[0];
231 if (x <= 0.0) return pdfGradient;
232 const Scalar pdf = computePDF(point);
233 pdfGradient[0] = -(std::log(2.0) + std::log(x) + SpecFunc::DiGamma(0.5 * nu_)) * pdf;
234 pdfGradient[1] = 0.5 * (0.5 / x - nu_) * pdf;
235 return pdfGradient;
236 }
237
238 /* Get the CDFGradient of the distribution */
computeCDFGradient(const Point & point) const239 Point InverseChiSquare::computeCDFGradient(const Point & point) const
240 {
241 if (point.getDimension() != 1) throw InvalidArgumentException(HERE) << "Error: the given point must have dimension=1, here dimension=" << point.getDimension();
242
243 Point cdfGradient(2, 0.0);
244 const Scalar x = point[0];
245 if (x <= 0.0) return cdfGradient;
246 const Scalar lambdaXInverse = 0.5 / x;
247 const Scalar pdf = computePDF(x);
248 const Scalar eps = std::pow(cdfEpsilon_, 1.0 / 3.0);
249 cdfGradient[0] = (DistFunc::pGamma(0.5 * nu_ + eps, lambdaXInverse, true) - DistFunc::pGamma(0.5 * nu_ - eps, lambdaXInverse, true)) / (2.0 * eps);
250 cdfGradient[1] = 0.5 * pdf * x;
251 return cdfGradient;
252 }
253
254 /* Get the quantile of the distribution */
computeScalarQuantile(const Scalar prob,const Bool tail) const255 Scalar InverseChiSquare::computeScalarQuantile(const Scalar prob,
256 const Bool tail) const
257 {
258 return 0.5 / DistFunc::qGamma(0.5 * nu_, prob, !tail);
259 }
260
261 /* Compute the mean of the distribution */
computeMean() const262 void InverseChiSquare::computeMean() const
263 {
264 if (!(nu_ > 2.0)) throw NotDefinedException(HERE) << "InverseChiSquare mean is defined only for nu > 2, here nu=" << nu_;
265 mean_ = Point(1, 1.0 / (nu_ - 2.0));
266 isAlreadyComputedMean_ = true;
267 }
268
269 /* Get the standard deviation of the distribution */
getStandardDeviation() const270 Point InverseChiSquare::getStandardDeviation() const
271 {
272 if (!(nu_ > 4.0)) throw NotDefinedException(HERE) << "InverseChiSquare standard deviation is defined only for nu > 4, here nu=" << nu_;
273 return Point(1, std::sqrt(getCovariance()(0, 0)));
274 }
275
276 /* Get the skewness of the distribution */
getSkewness() const277 Point InverseChiSquare::getSkewness() const
278 {
279 if (!(nu_ > 6.0)) throw NotDefinedException(HERE) << "InverseChiSquare skewness is defined only for nu > 6, here nu=" << nu_;
280 return Point(1, 8.0 * std::sqrt(0.5 * nu_ - 2.0) / (nu_ - 6.0));
281 }
282
283 /* Get the kurtosis of the distribution */
getKurtosis() const284 Point InverseChiSquare::getKurtosis() const
285 {
286 if (!(nu_ > 8.0)) throw NotDefinedException(HERE) << "InverseChiSquare kurtosis is defined only for nu > 8, here nu=" << nu_;
287 return Point(1, 12.0 * (0.5 * nu_ * (0.5 * nu_ + 3.0) - 10.0) / ((nu_ - 6.0) * (nu_ - 8.0)));
288 }
289
290 /* Get the moments of the standardized distribution */
getStandardMoment(const UnsignedInteger n) const291 Point InverseChiSquare::getStandardMoment(const UnsignedInteger n) const
292 {
293 if (nu_ <= 2.0 * n) throw NotDefinedException(HERE) << "InverseChiSquare standard moment of order " << 2.0 * n << " is defined only for nu > " << 2.0 * n << ", here k=" << nu_;
294 return Point(1, std::exp(SpecFunc::LogGamma(0.5 * nu_ - n) - SpecFunc::LogGamma(0.5 * nu_)));
295 }
296
297 /* Get the standard representative in the parametric family, associated with the standard moments */
getStandardRepresentative() const298 Distribution InverseChiSquare::getStandardRepresentative() const
299 {
300 return new InverseChiSquare(nu_);
301 }
302
303 /* Compute the covariance of the distribution */
computeCovariance() const304 void InverseChiSquare::computeCovariance() const
305 {
306 if (!(nu_ > 4.0)) throw NotDefinedException(HERE) << "InverseChiSquare covariance is defined only for nu > 4, here nu=" << nu_;
307 covariance_ = CovarianceMatrix(1);
308 covariance_(0, 0) = 2.0 / ((nu_ - 2.0) * (nu_ - 2.0) * (nu_ - 4.0));
309 isAlreadyComputedCovariance_ = true;
310 }
311
312 /* Parameters value accessor */
getParameter() const313 Point InverseChiSquare::getParameter() const
314 {
315 return Point(1, nu_);
316 }
317
setParameter(const Point & parameter)318 void InverseChiSquare::setParameter(const Point & parameter)
319 {
320 if (parameter.getSize() != 1) throw InvalidArgumentException(HERE) << "Error: expected 1 value, got " << parameter.getSize();
321 const Scalar w = getWeight();
322 *this = InverseChiSquare(parameter[0]);
323 setWeight(w);
324 }
325
326 /* Parameters description accessor */
getParameterDescription() const327 Description InverseChiSquare::getParameterDescription() const
328 {
329 return Description(1, "nu");
330 }
331
332 /* Method save() stores the object through the StorageManager */
save(Advocate & adv) const333 void InverseChiSquare::save(Advocate & adv) const
334 {
335 ContinuousDistribution::save(adv);
336 adv.saveAttribute( "nu_", nu_ );
337 adv.saveAttribute( "normalizationFactor_", normalizationFactor_ );
338 }
339
340 /* Method load() reloads the object from the StorageManager */
load(Advocate & adv)341 void InverseChiSquare::load(Advocate & adv)
342 {
343 ContinuousDistribution::load(adv);
344 adv.loadAttribute( "nu_", nu_ );
345 adv.loadAttribute( "normalizationFactor_", normalizationFactor_ );
346 computeRange();
347 }
348
349 END_NAMESPACE_OPENTURNS
350