20
20
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
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 $
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)
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
156
157
NumericalScalar x(point[0]);
157
158
if (x <= a_) return (tail ? 1.0 : 0.0);
202
203
/* Get the quantile of the distribution */
203
204
NumericalScalar Beta::computeScalarQuantile(const NumericalScalar prob,
204
const NumericalScalar initialGuess,
205
const NumericalScalar initialStep,
206
206
const NumericalScalar precision) const
208
return a_ + (b_ - a_) * DistFunc::qBeta(r_, t_ - r_, prob);
208
return a_ + (b_ - a_) * DistFunc::qBeta(r_, t_ - r_, prob, tail);
211
211
/* Get the roughness, i.e. the L2-norm of the PDF */
218
218
/* Get the mean of the distribution */
219
Beta::NumericalPoint Beta::getMean() const throw(NotDefinedException)
219
Beta::NumericalPoint Beta::getMean() const /* throw(NotDefinedException) */
221
221
return NumericalPoint(1, a_ + (b_ - a_) * r_ / t_);
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) */
227
227
return NumericalPoint(1, getSigma());
230
230
/* Get the skewness of the distribution */
231
Beta::NumericalPoint Beta::getSkewness() const throw(NotDefinedException)
231
Beta::NumericalPoint Beta::getSkewness() const /* throw(NotDefinedException) */
233
233
return NumericalPoint(1, 2.0 * (t_ - 2.0 * r_) / (t_ + 2.0) * sqrt((t_ + 1.0) / (r_ * (t_ - r_))));
236
236
/* Get the kurtosis of the distribution */
237
Beta::NumericalPoint Beta::getKurtosis() const throw(NotDefinedException)
237
Beta::NumericalPoint Beta::getKurtosis() const /* throw(NotDefinedException) */
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_)));
242
242
/* Get the covariance of the distribution */
243
Beta::CovarianceMatrix Beta::getCovariance() const throw(NotDefinedException)
243
Beta::CovarianceMatrix Beta::getCovariance() const /* throw(NotDefinedException) */
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;
251
/* Get the moments of the standardized distribution */
252
Beta::NumericalPoint Beta::getStandardMoment(const UnsignedLong n) const
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)
265
term *= -2.0 * (r_ + j) / (t_ + j) * (n - j) / (j + 1);
268
return NumericalPoint(1, value);
251
271
/* Parameters value and description accessor */
252
272
Beta::NumericalPointWithDescriptionCollection Beta::getParametersCollection() const
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
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_ );
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)
387
407
NonEllipticalDistribution::load(adv);
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;
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_ );