1 // -*- C++ -*-
2 /**
3 * @brief The ChiSquare distribution, ie the Gamma(nu/2,1/2) 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/ChiSquare.hxx"
23 #include "openturns/Gamma.hxx"
24 #include "openturns/RandomGenerator.hxx"
25 #include "openturns/SpecFunc.hxx"
26 #include "openturns/DistFunc.hxx"
27 #include "openturns/PersistentObjectFactory.hxx"
28
29 BEGIN_NAMESPACE_OPENTURNS
30
31 CLASSNAMEINIT(ChiSquare)
32
33 static const Factory<ChiSquare> Factory_ChiSquare;
34
35 /* Default constructor */
ChiSquare()36 ChiSquare::ChiSquare()
37 : ContinuousDistribution()
38 , nu_(0.0)
39 , normalizationFactor_(0.0)
40 {
41 setName("ChiSquare");
42 // This call also call update() and computeRange()
43 setNu(1.0);
44 setDimension(1);
45 computeRange();
46 }
47
48 /* Parameters constructor */
ChiSquare(const Scalar nu)49 ChiSquare::ChiSquare(const Scalar nu)
50 : ContinuousDistribution()
51 , nu_(0.0)
52 , normalizationFactor_(0.0)
53 {
54 setName("ChiSquare");
55 // This call set also call computeRange() and update().
56 setNu(nu);
57 setDimension(1);
58 }
59
60 /* Comparison operator */
operator ==(const ChiSquare & other) const61 Bool ChiSquare::operator ==(const ChiSquare & other) const
62 {
63 if (this == &other) return true;
64 return nu_ == other.nu_;
65 }
66
equals(const DistributionImplementation & other) const67 Bool ChiSquare::equals(const DistributionImplementation & other) const
68 {
69 const ChiSquare* p_other = dynamic_cast<const ChiSquare*>(&other);
70 return p_other && (*this == *p_other);
71 }
72
73 /* String converter */
__repr__() const74 String ChiSquare::__repr__() const
75 {
76 OSS oss(true);
77 oss << "class=" << ChiSquare::GetClassName()
78 << " name=" << getName()
79 << " dimension=" << getDimension()
80 << " nu=" << nu_;
81 return oss;
82 }
83
__str__(const String &) const84 String ChiSquare::__str__(const String & ) const
85 {
86 OSS oss(false);
87 oss << getClassName() << "(nu = " << nu_ << ")";
88 return oss;
89 }
90
91 /* Nu accessor */
setNu(const Scalar nu)92 void ChiSquare::setNu(const Scalar nu)
93 {
94 if (!(nu > 0.0)) throw InvalidArgumentException(HERE) << "Nu MUST be positive";
95 if (nu != nu_)
96 {
97 nu_ = nu;
98 computeRange();
99 update();
100 }
101 }
102
getNu() const103 Scalar ChiSquare::getNu() const
104 {
105 return nu_;
106 }
107
108
109 /* Virtual constructor */
clone() const110 ChiSquare * ChiSquare::clone() const
111 {
112 return new ChiSquare(*this);
113 }
114
115 /* Compute the numerical range of the distribution given the parameters values */
computeRange()116 void ChiSquare::computeRange()
117 {
118 setRange(Gamma(0.5, 0.5 * nu_, 0.0).getRange());
119 }
120
121 /** Update the derivative attributes */
update()122 void ChiSquare::update()
123 {
124 normalizationFactor_ = -0.5 * nu_ * M_LN2 - SpecFunc::LnGamma(0.5 * nu_);
125 isAlreadyComputedMean_ = false;
126 isAlreadyComputedCovariance_ = false;
127 }
128
129
130 /* Get one realization of the distribution */
getRealization() const131 Point ChiSquare::getRealization() const
132 {
133 return Point(1, 2.0 * DistFunc::rGamma(0.5 * nu_));
134 }
135
136
137 /* Get the DDF of the distribution */
computeDDF(const Point & point) const138 Point ChiSquare::computeDDF(const Point & point) const
139 {
140 if (point.getDimension() != 1) throw InvalidArgumentException(HERE) << "Error: the given point must have dimension=1, here dimension=" << point.getDimension();
141
142 const Scalar x = point[0];
143 if (x <= 0.0) return Point(1, 0.0);
144 return Point(1, ((0.5 * nu_ - 1.0) / x - 0.5) * computePDF(point));
145 }
146
147
148 /* Get the PDF of the distribution */
computePDF(const Point & point) const149 Scalar ChiSquare::computePDF(const Point & point) const
150 {
151 if (point.getDimension() != 1) throw InvalidArgumentException(HERE) << "Error: the given point must have dimension=1, here dimension=" << point.getDimension();
152
153 const Scalar x = point[0];
154 if (x <= 0.0) return 0.0;
155 return std::exp(computeLogPDF(point));
156 }
157
computeLogPDF(const Point & point) const158 Scalar ChiSquare::computeLogPDF(const Point & point) const
159 {
160 if (point.getDimension() != 1) throw InvalidArgumentException(HERE) << "Error: the given point must have dimension=1, here dimension=" << point.getDimension();
161
162 const Scalar x = point[0];
163 if (x <= 0.0) return SpecFunc::LowestScalar;
164 return normalizationFactor_ + (0.5 * nu_ - 1) * std::log(x) - 0.5 * x;
165 }
166
167 /* Get the CDF of the distribution */
computeCDF(const Point & point) const168 Scalar ChiSquare::computeCDF(const Point & point) const
169 {
170 if (point.getDimension() != 1) throw InvalidArgumentException(HERE) << "Error: the given point must have dimension=1, here dimension=" << point.getDimension();
171
172 const Scalar x = point[0];
173 // No test here as the CDF is continuous for all nu_
174 if (x <= 0.0) return 0.0;
175 return DistFunc::pGamma(0.5 * nu_, 0.5 * x);
176 }
177
computeComplementaryCDF(const Point & point) const178 Scalar ChiSquare::computeComplementaryCDF(const Point & point) const
179 {
180 if (point.getDimension() != 1) throw InvalidArgumentException(HERE) << "Error: the given point must have dimension=1, here dimension=" << point.getDimension();
181
182 const Scalar x = point[0];
183 // No test here as the CDF is continuous for all nu_
184 if (x <= 0.0) return 1.0;
185 return DistFunc::pGamma(0.5 * nu_, 0.5 * x, true);
186 }
187
188 /* Compute the entropy of the distribution */
computeEntropy() const189 Scalar ChiSquare::computeEntropy() const
190 {
191 return 0.5 * nu_ + M_LN2 + SpecFunc::LogGamma(0.5 * nu_) + (1.0 - 0.5 * nu_) * SpecFunc::Psi(0.5 * nu_);
192 }
193
194 /* Get the characteristic function of the distribution, i.e. phi(u) = E(exp(I*u*X)) */
computeCharacteristicFunction(const Scalar x) const195 Complex ChiSquare::computeCharacteristicFunction(const Scalar x) const
196 {
197 return std::pow(Complex(1.0, -2.0 * x), -0.5 * nu_);
198 }
199
computeLogCharacteristicFunction(const Scalar x) const200 Complex ChiSquare::computeLogCharacteristicFunction(const Scalar x) const
201 {
202 return -0.5 * nu_ * std::log(Complex(1.0, -2.0 * x));
203 }
204
205 /* Get the PDFGradient of the distribution */
computePDFGradient(const Point & point) const206 Point ChiSquare::computePDFGradient(const Point & point) const
207 {
208 if (point.getDimension() != 1) throw InvalidArgumentException(HERE) << "Error: the given point must have dimension=1, here dimension=" << point.getDimension();
209
210 Point pdfGradient(1, 0.0);
211 const Scalar x = point[0];
212 if (x <= 0.0) return pdfGradient;
213 Scalar pdf = computePDF(point);
214 pdfGradient[0] = 0.5 * (std::log(0.5 * x) - SpecFunc::Psi(0.5 * nu_)) * pdf;
215 return pdfGradient;
216 }
217
218 /* Get the CDFGradient of the distribution */
computeCDFGradient(const Point & point) const219 Point ChiSquare::computeCDFGradient(const Point & point) const
220 {
221 if (point.getDimension() != 1) throw InvalidArgumentException(HERE) << "Error: the given point must have dimension=1, here dimension=" << point.getDimension();
222
223 Point cdfGradient(1, 0.0);
224 const Scalar x = point[0];
225 if (x <= 0.0) return cdfGradient;
226 Scalar eps = std::pow(cdfEpsilon_, 1.0 / 3.0);
227 cdfGradient[0] = (DistFunc::pGamma(0.5 * (nu_ + eps), 0.5 * x) - DistFunc::pGamma(0.5 * (nu_ - eps), 0.5 * x)) / (2.0 * eps);
228 return cdfGradient;
229 }
230
231 /* Get the quantile of the distribution */
computeScalarQuantile(const Scalar prob,const Bool tail) const232 Scalar ChiSquare::computeScalarQuantile(const Scalar prob,
233 const Bool tail) const
234 {
235 return 2.0 * DistFunc::qGamma(0.5 * nu_, prob, tail);
236 }
237
238 /* Compute the mean of the distribution */
computeMean() const239 void ChiSquare::computeMean() const
240 {
241 mean_ = Point(1, nu_);
242 isAlreadyComputedMean_ = true;
243 }
244
245 /* Get the standard deviation of the distribution */
getStandardDeviation() const246 Point ChiSquare::getStandardDeviation() const
247 {
248 return Point(1, std::sqrt(2.0 * nu_));
249 }
250
251 /* Get the skewness of the distribution */
getSkewness() const252 Point ChiSquare::getSkewness() const
253 {
254 return Point(1, std::sqrt(8.0 / nu_));
255 }
256
257 /* Get the kurtosis of the distribution */
getKurtosis() const258 Point ChiSquare::getKurtosis() const
259 {
260 return Point(1, 3.0 + 12.0 / nu_);
261 }
262
263 /* Compute the covariance of the distribution */
computeCovariance() const264 void ChiSquare::computeCovariance() const
265 {
266 covariance_ = CovarianceMatrix(1);
267 covariance_(0, 0) = 2.0 * nu_;
268 isAlreadyComputedCovariance_ = true;
269 }
270
271 /* Get the moments of the standardized distribution */
getStandardMoment(const UnsignedInteger n) const272 Point ChiSquare::getStandardMoment(const UnsignedInteger n) const
273 {
274 return Point(1, std::exp(SpecFunc::LnGamma(n + 0.5 * nu_) - SpecFunc::LnGamma(0.5 * nu_)));
275 }
276
277 /* Get the standard representative in the parametric family, associated with the standard moments */
getStandardRepresentative() const278 Distribution ChiSquare::getStandardRepresentative() const
279 {
280 return Gamma(0.5 * nu_, 1.0, 0.0);
281 }
282
283 /* Parameters value accessor */
getParameter() const284 Point ChiSquare::getParameter() const
285 {
286 return Point(1, nu_);
287 }
288
setParameter(const Point & parameter)289 void ChiSquare::setParameter(const Point & parameter)
290 {
291 if (parameter.getSize() != 1) throw InvalidArgumentException(HERE) << "Error: expected 1 value, got " << parameter.getSize();
292 const Scalar w = getWeight();
293 *this = ChiSquare(parameter[0]);
294 setWeight(w);
295 }
296
297 /* Parameters description accessor */
getParameterDescription() const298 Description ChiSquare::getParameterDescription() const
299 {
300 return Description(1, "nu");
301 }
302
303 /* Method save() stores the object through the StorageManager */
save(Advocate & adv) const304 void ChiSquare::save(Advocate & adv) const
305 {
306 ContinuousDistribution::save(adv);
307 adv.saveAttribute( "nu_", nu_ );
308 adv.saveAttribute( "normalizationFactor_", normalizationFactor_ );
309 }
310
311 /* Method load() reloads the object from the StorageManager */
load(Advocate & adv)312 void ChiSquare::load(Advocate & adv)
313 {
314 ContinuousDistribution::load(adv);
315 adv.loadAttribute( "nu_", nu_ );
316 adv.loadAttribute( "normalizationFactor_", normalizationFactor_ );
317 computeRange();
318 }
319
320
321
322
323 END_NAMESPACE_OPENTURNS
324