~ubuntu-branches/ubuntu/maverick/openturns/maverick

« back to all changes in this revision

Viewing changes to lib/src/Uncertainty/Distribution/Beta.cxx

  • Committer: Bazaar Package Importer
  • Author(s): Fabrice Coutadeur
  • Date: 2010-05-10 17:27:55 UTC
  • mfrom: (1.1.4 upstream) (5.1.5 sid)
  • Revision ID: james.westby@ubuntu.com-20100510172755-cb5ynskknqqi5rhp
Tags: 0.13.2-2ubuntu1
* Merge with Debian testing. No changes left.
* ubuntu_fix-python-2.6.patch: fix detection of python 2.6 libs, to not use
  LOCALMODLIBS. This pulls a dependency on SSL and makes the package FTBFS.

Show diffs side-by-side

added added

removed removed

Lines of Context:
3
3
 *  @file  Beta.cxx
4
4
 *  @brief The Beta distribution
5
5
 *
6
 
 *  (C) Copyright 2005-2007 EDF-EADS-Phimeca
 
6
 *  (C) Copyright 2005-2010 EDF-EADS-Phimeca
7
7
 *
8
8
 *  This library is free software; you can redistribute it and/or
9
9
 *  modify it under the terms of the GNU Lesser General Public
20
20
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
21
21
 *
22
22
 *  @author: $LastChangedBy: dutka $
23
 
 *  @date:   $LastChangedDate: 2009-05-28 14:47:53 +0200 (jeu. 28 mai 2009) $
24
 
 *  Id:      $Id: Beta.cxx 1262 2009-05-28 12:47:53Z dutka $
 
23
 *  @date:   $LastChangedDate: 2010-02-04 16:44:49 +0100 (jeu. 04 févr. 2010) $
 
24
 *  Id:      $Id: Beta.cxx 1473 2010-02-04 15:44:49Z dutka $
25
25
 */
26
26
#include <cmath>
27
27
#include "Beta.hxx"
57
57
                 const NumericalScalar a,
58
58
                 const NumericalScalar b,
59
59
                 const ParameterSet set)
60
 
        throw (InvalidArgumentException)
 
60
        /* throw (InvalidArgumentException) */
61
61
        : NonEllipticalDistribution("Beta"),
62
62
          r_(0.), t_(0.), a_(a), b_(b), normalizationFactor_(0.0)
63
63
      {
151
151
 
152
152
 
153
153
      /* Get the CDF of the distribution */
154
 
      NumericalScalar Beta::computeCDF(const NumericalPoint & point, const Bool tail) const
 
154
      NumericalScalar Beta::computeCDF(const NumericalPoint & point,
 
155
                                       const Bool tail) const
155
156
      {
156
157
        NumericalScalar x(point[0]);
157
158
        if (x <= a_) return (tail ? 1.0 : 0.0);
201
202
 
202
203
      /* Get the quantile of the distribution */
203
204
      NumericalScalar Beta::computeScalarQuantile(const NumericalScalar prob,
204
 
                                                  const NumericalScalar initialGuess,
205
 
                                                  const NumericalScalar initialStep,
 
205
                                                  const Bool tail,
206
206
                                                  const NumericalScalar precision) const
207
207
      {
208
 
        return a_ + (b_ - a_) * DistFunc::qBeta(r_, t_ - r_, prob);
 
208
        return a_ + (b_ - a_) * DistFunc::qBeta(r_, t_ - r_, prob, tail);
209
209
      }
210
210
 
211
211
      /* Get the roughness, i.e. the L2-norm of the PDF */
216
216
      }
217
217
 
218
218
      /* Get the mean of the distribution */
219
 
      Beta::NumericalPoint Beta::getMean() const throw(NotDefinedException)
 
219
      Beta::NumericalPoint Beta::getMean() const /* throw(NotDefinedException) */
220
220
      {
221
221
        return NumericalPoint(1, a_ + (b_ - a_) * r_ / t_);
222
222
      }
223
223
 
224
224
      /* Get the standard deviation of the distribution */
225
 
      Beta::NumericalPoint Beta::getStandardDeviation() const throw(NotDefinedException)
 
225
      Beta::NumericalPoint Beta::getStandardDeviation() const /* throw(NotDefinedException) */
226
226
      {
227
227
        return NumericalPoint(1, getSigma());
228
228
      }
229
229
 
230
230
      /* Get the skewness of the distribution */
231
 
      Beta::NumericalPoint Beta::getSkewness() const throw(NotDefinedException)
 
231
      Beta::NumericalPoint Beta::getSkewness() const /* throw(NotDefinedException) */
232
232
      {
233
233
        return NumericalPoint(1, 2.0 * (t_ - 2.0 * r_) / (t_ + 2.0) * sqrt((t_ + 1.0) / (r_ * (t_ - r_))));
234
234
      }
235
235
 
236
236
      /* Get the kurtosis of the distribution */
237
 
      Beta::NumericalPoint Beta::getKurtosis() const throw(NotDefinedException)
 
237
      Beta::NumericalPoint Beta::getKurtosis() const /* throw(NotDefinedException) */
238
238
      {
239
239
        return NumericalPoint(1, 3.0 * (1.0 + t_) * (2.0 * t_ * t_ + r_ * (t_ - 6.0) * (t_ - r_))/ (r_ * (t_ - r_) * (3.0 + t_) * (2.0 + t_)));
240
240
      }
241
241
 
242
242
      /* Get the covariance of the distribution */
243
 
      Beta::CovarianceMatrix Beta::getCovariance() const throw(NotDefinedException)
 
243
      Beta::CovarianceMatrix Beta::getCovariance() const /* throw(NotDefinedException) */
244
244
      {
245
245
        CovarianceMatrix covariance(1);
246
 
        NumericalScalar eta((b_ - a_) / t_);
 
246
        const NumericalScalar eta((b_ - a_) / t_);
247
247
        covariance(0, 0) = eta * eta * r_ * (t_ - r_) / (t_ + 1.0);
248
248
        return covariance;
249
249
      }
250
250
 
 
251
      /* Get the moments of the standardized distribution */
 
252
      Beta::NumericalPoint Beta::getStandardMoment(const UnsignedLong n) const
 
253
      {
 
254
        if (n == 0) return NumericalPoint(1, 1.0);
 
255
        // The raw moments of the Beta distribution have a closed form in
 
256
        // term of the Gamma function for a=0, b=1:
 
257
        // m(k) = Gamma(t) * Gamma(r + n) / (Gamma(r) * Gamma(t + n))
 
258
        // Here, the standard distribution corresponds to a=-1, b=1 so
 
259
        // we have to perform a little algebra to get the needed values
 
260
        NumericalScalar term(1.0);
 
261
        if (n % 2 == 1) term = -term;
 
262
        NumericalScalar value(term);
 
263
        for (UnsignedLong j = 1; j < n; ++j)
 
264
          {
 
265
            term *= -2.0 * (r_ + j) / (t_ + j) * (n - j) / (j + 1);
 
266
            value += term;
 
267
          }
 
268
        return NumericalPoint(1, value);
 
269
      }
 
270
 
251
271
      /* Parameters value and description accessor */
252
272
      Beta::NumericalPointWithDescriptionCollection Beta::getParametersCollection() const
253
273
      {
283
303
 
284
304
      /* R accessor */
285
305
      void Beta::setR(const NumericalScalar r)
286
 
        throw(InvalidArgumentException)
 
306
        /* throw(InvalidArgumentException) */
287
307
      {
288
308
        if (r <= 0.) throw InvalidArgumentException(HERE) << "R MUST be positive";
289
309
        r_ = r;
298
318
 
299
319
      /* T accessor */
300
320
      void Beta::setT(const NumericalScalar t)
301
 
        throw(InvalidArgumentException)
 
321
        /* throw(InvalidArgumentException) */
302
322
      {
303
323
        if (t <= r_) throw InvalidArgumentException(HERE) << "T MUST be greater than r, here t=" << t << " and r=" << r_;
304
324
        t_ = t;
314
334
      /* RT accessor */
315
335
      void Beta::setRT(const NumericalScalar r,
316
336
                       const NumericalScalar t)
317
 
        throw(InvalidArgumentException)
 
337
        /* throw(InvalidArgumentException) */
318
338
      {
319
339
        if (r <= 0.) throw InvalidArgumentException(HERE) << "R MUST be positive";
320
340
        r_ = r;
371
391
      }
372
392
 
373
393
      /* Method save() stores the object through the StorageManager */
374
 
      void Beta::save(const StorageManager::Advocate & adv) const
 
394
      void Beta::save(StorageManager::Advocate & adv) const
375
395
      {
376
396
        NonEllipticalDistribution::save(adv);
377
 
        adv.writeValue("r_", r_);
378
 
        adv.writeValue("t_", t_);
379
 
        adv.writeValue("a_", a_);
380
 
        adv.writeValue("b_", b_);
381
 
        adv.writeValue("normalizationFactor_", normalizationFactor_);
 
397
        adv.saveAttribute( "r_", r_ );
 
398
        adv.saveAttribute( "t_", t_ );
 
399
        adv.saveAttribute( "a_", a_ );
 
400
        adv.saveAttribute( "b_", b_ );
 
401
        adv.saveAttribute( "normalizationFactor_", normalizationFactor_ );
382
402
      }
383
403
 
384
404
      /* Method load() reloads the object from the StorageManager */
385
 
      void Beta::load(const StorageManager::Advocate & adv)
 
405
      void Beta::load(StorageManager::Advocate & adv)
386
406
      {
387
407
        NonEllipticalDistribution::load(adv);
388
 
 
389
 
        String name;
390
 
        NumericalScalar value;
391
 
        StorageManager::List objList = adv.getList(StorageManager::NumericalScalarEntity);
392
 
        for(objList.firstValueToRead(); objList.moreValuesToRead(); objList.nextValueToRead()) {
393
 
          if (objList.readValue(name, value)) {
394
 
            if (name == "r_") r_ = value;
395
 
            if (name == "t_") t_ = value;
396
 
            if (name == "a_") a_ = value;
397
 
            if (name == "b_") b_ = value;
398
 
            if (name == "normalizationFactor_") normalizationFactor_ = value;
399
 
          }
400
 
        }
 
408
        adv.loadAttribute( "r_", r_ );
 
409
        adv.loadAttribute( "t_", t_ );
 
410
        adv.loadAttribute( "a_", a_ );
 
411
        adv.loadAttribute( "b_", b_ );
 
412
        adv.loadAttribute( "normalizationFactor_", normalizationFactor_ );
401
413
        computeRange();
402
414
      }
403
415