Page MenuHomeHEPForge
Files F10881666

No OneTemporary
Actions

View Options

Index: trunk/npstat/stat/SeriesCGF1D.hh
===================================================================
--- trunk/npstat/stat/SeriesCGF1D.hh (revision 754)
+++ trunk/npstat/stat/SeriesCGF1D.hh (revision 755)
@@ -1,59 +1,65 @@
#ifndef NPSTAT_SERIESCGF1D_HH_
#define NPSTAT_SERIESCGF1D_HH_
/*!
// \file SeriesCGF1D.hh
//
// \brief Univariate cumulant generating function defined by cumulants
//
// Author: I. Volobouev
//
// December 2019
*/
#include <vector>
#include "npstat/stat/AbsCGF1D.hh"
namespace npstat {
class SeriesCGF1D : public AbsCGF1D
{
public:
/**
// At least the first two cumulants (mean and variance)
// must be provided in the "cumulants" array.
// cumulants[0] must contain the mean, cumulants[1] the
// variance, etc. "nCumulants" specifies the length of
// the array (and also the highest cumulant order).
// The variance must be positive (i.e., the degenerate
// distribution of a single point mass is excluded from
// consideration). The parameters "smin" and "smax"
// should be chosen in such a way that the CGF remains
// convex on the (smin, smax) interval.
*/
SeriesCGF1D(const double* cumulants, unsigned nCumulants,
double smin, double smax);
inline virtual SeriesCGF1D* clone() const
{return new SeriesCGF1D(*this);}
inline virtual ~SeriesCGF1D() {}
virtual SeriesCGF1D* shiftAndScale(double mu, double sigma) const;
inline virtual double smin() const {return smin_;}
inline virtual double smax() const {return smax_;}
virtual double derivative(unsigned order, double s) const;
+ /**
+ // Check whether the CGF is convex on the interval (a, b].
+ // This method ignores smin and smax.
+ */
+ virtual bool isConvex(double a, double b) const;
+
protected:
virtual bool isEqual(const AbsCGF1D& r) const;
private:
std::vector<double> cumulants_;
double smin_;
double smax_;
};
}
#endif // NPSTAT_SERIESCGF1D_HH_
Index: trunk/npstat/stat/SeriesCGF1D.cc
===================================================================
--- trunk/npstat/stat/SeriesCGF1D.cc (revision 754)
+++ trunk/npstat/stat/SeriesCGF1D.cc (revision 755)
@@ -1,92 +1,108 @@
#include <cmath>
#include <cassert>
#include <stdexcept>
#include "npstat/stat/SeriesCGF1D.hh"
#include "npstat/rng/permutation.hh"
+#include "npstat/nm/Poly1D.hh"
namespace npstat {
SeriesCGF1D::SeriesCGF1D(const double* cumulants,
const unsigned nCumulants,
const double i_smin, const double i_smax)
: cumulants_(cumulants, cumulants+nCumulants),
smin_(i_smin),
smax_(i_smax)
{
assert(cumulants);
if (nCumulants < 2U) throw std::invalid_argument(
- "In npstat::SeriesCGF1D constructor : not enough cumulants");
+ "In npstat::SeriesCGF1D constructor: not enough cumulants");
if (cumulants[1] <= 0.0) throw std::invalid_argument(
- "In npstat::SeriesCGF1D constructor : variance must be positive");
+ "In npstat::SeriesCGF1D constructor: variance must be positive");
if (!(smin_ < 0.0)) throw std::invalid_argument(
- "In npstat::SeriesCGF1D constructor : invalid support infimum");
+ "In npstat::SeriesCGF1D constructor: invalid support infimum");
if (!(smax_ > 0.0)) throw std::invalid_argument(
- "In npstat::SeriesCGF1D constructor : invalid support supremum");
+ "In npstat::SeriesCGF1D constructor: invalid support supremum");
}
bool SeriesCGF1D::isEqual(const AbsCGF1D& other) const
{
const SeriesCGF1D& r = static_cast<const SeriesCGF1D&>(other);
return cumulants_ == r.cumulants_ && smin_ == r.smin_ && smax_ == r.smax_;
}
double SeriesCGF1D::derivative(const unsigned order, const double s) const
{
if (s == 0.0)
{
if (order)
{
const unsigned idx = order - 1U;
if (idx < cumulants_.size())
return cumulants_[idx];
}
}
else if (!(smin_ < s && s < smax_)) throw std::invalid_argument(
- "In npstat::SeriesCGF1D::derivative : argument out of range");
+ "In npstat::SeriesCGF1D::derivative: argument out of range");
else
{
const unsigned degree = cumulants_.size();
if (order <= degree)
{
const double* a = &cumulants_[0];
long double sum = 0.0L;
if (order)
{
const int shift = order-1U;
for (int deg=degree-order; deg>=0; --deg)
{
sum *= s;
sum += a[deg+shift]/ldfactorial(deg);
}
}
else
{
for (unsigned deg=degree; deg>0U; --deg)
{
sum += a[deg-1]/ldfactorial(deg);
sum *= s;
}
}
return sum;
}
}
return 0.0;
}
SeriesCGF1D* SeriesCGF1D::shiftAndScale(const double mu,
const double sigma) const
{
if (!(sigma > 0.0)) throw std::invalid_argument(
- "In npstat::SeriesCGF1D::shiftAndScale : scale must be positive");
+ "In npstat::SeriesCGF1D::shiftAndScale: scale must be positive");
std::vector<double> cum(cumulants_);
const unsigned sz = cum.size();
cum[0] -= mu;
if (sigma != 1.0)
for (unsigned i=0; i<sz; ++i)
cum[i] /= pow(sigma, i+1U);
// Rethink this. It is not obvious how the
// support should be modified.
return new SeriesCGF1D(&cum[0], sz, smin_, smax_);
}
+
+ bool SeriesCGF1D::isConvex(const double a, const double b) const
+ {
+ if (a >= b) throw std::invalid_argument(
+ "In npstat::SeriesCGF1D::sisConvex: invalid interval");
+ const unsigned degree = cumulants_.size();
+ std::vector<long double> coeffs(degree + 1U);
+ coeffs[0] = 0.0L;
+ for (unsigned deg=0; deg<degree; ++deg)
+ coeffs[deg+1] = cumulants_[deg]/ldfactorial(deg+1U);
+ const Poly1D cgf0(coeffs);
+ const Poly1D& cgf1(cgf0.derivative());
+ const Poly1D& cgf2(cgf1.derivative());
+ return cgf2((a + b)/2.0) > 0.0L && !cgf2.nRoots(a, b);
+ }
}

File Metadata

Mime Type
text/x-diff
Expires
Sat, May 3, 6:40 AM (22 h, 50 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4983113
Default Alt Text
(6 KB)

Event Timeline