EllipticalDistributions.hh

Go to the documentation of this file.

#ifndef NPSTAT_ELLIPTICALDISTRIBUTIONS_HH_
#define NPSTAT_ELLIPTICALDISTRIBUTIONS_HH_
 
/*!
// \file EllipticalDistributions.hh
//
// \brief A few concrete examples of multivariate elliptical ditributions
//        with simple parameters
//
// Author: I. Volobouev
//
// September 2022
*/
 
#include "npstat/stat/EllipticalDistribution.hh"
#include "npstat/stat/Distributions1D.hh"
#include "npstat/stat/PowerRatio1D.hh"
#include "npstat/stat/PowerTransform1D.hh"
#include "npstat/stat/TransformedDistribution1D.hh"
 
// The distributions below are taken from Chapter 3 of the
// monograph "Symmetric Multivariate and Related Distributions"
// by K.-T. Fang, S. Kotz and K. W. Ng, CRC Press, 1990.
namespace npstat {
    /**
    // EllipticalNormal is the same distribution as GaussND.
    // The transformation matrix in the constructor is just
    // the square root of the desired covariance matrix.
    */
    class EllipticalNormal : public EllipticalDistribution
    {
    public:
        inline EllipticalNormal(const double* location1, unsigned dim1,
                                const Matrix<double>& transformationMatrix)
            : EllipticalDistribution(location1, dim1,
                                     transformationMatrix,
                                     Exponential1D(0.0, 2.0),
                                     Gamma1D(0.0, 2.0, dim1/2.0)) {}
 
        inline virtual ~EllipticalNormal() {}
 
        // Methods needed for I/O
        inline virtual gs::ClassId classId() const {return gs::ClassId(*this);}
        inline virtual bool write(std::ostream& os) const
        {
            return EllipticalDistribution::classId().write(os) &&
                   EllipticalDistribution::write(os);
        }
 
        static inline const char* classname()
            {return "npstat::EllipticalNormal";}
        static inline unsigned version() {return 1;}
        static inline EllipticalNormal* read(
            const gs::ClassId& id, std::istream& in)
        {
            static const gs::ClassId current(
                gs::ClassId::makeId<EllipticalNormal>());
            current.ensureSameId(id);
            CPP11_auto_ptr<EllipticalDistribution> pb =
                gs::read_obj<EllipticalDistribution>(in);
            return new EllipticalNormal(*pb);
        }
 
    private:
        inline EllipticalNormal(const EllipticalDistribution& b)
            : EllipticalDistribution(b) {}
 
#ifdef SWIG
    public:
        inline EllipticalNormal(const double* location1, unsigned dim1,
                                const double* transform, unsigned nrows, unsigned ncols)
            : EllipticalDistribution(location1, dim1,
                                     transform, nrows, ncols,
                                     Exponential1D(0.0, 2.0),
                                     Gamma1D(0.0, 2.0, dim1/2.0)) {}
#endif // SWIG
    };
 
    class EllipticalKotz : public EllipticalDistribution
    {
    public:
        /** For this distribution we must have s > 0, r > 0, 2*N + dim1 > 2 */
        inline EllipticalKotz(const double* location1, unsigned dim1,
                              const Matrix<double>& transformationMatrix,
                              const unsigned N, const double r, const double s)
            : EllipticalDistribution(
                location1, dim1, transformationMatrix,
                TransformedDistribution1D(PowerTransform1D(r, s), Gamma1D(0.0, 1.0, N/s)),
                TransformedDistribution1D(PowerTransform1D(r, s), Gamma1D(0.0, 1.0, (N+dim1/2.0-1)/s))),
                N_(N), r_(r), s_(s) {calcNorm();}
 
        inline virtual ~EllipticalKotz() {}
 
        virtual double density(const double* location1, unsigned dim1) const;
 
        // Methods needed for I/O
        inline virtual gs::ClassId classId() const {return gs::ClassId(*this);}
        virtual bool write(std::ostream& os) const;
 
        static inline const char* classname()
            {return "npstat::EllipticalKotz";}
        static inline unsigned version() {return 1;}
        static EllipticalKotz* read(const gs::ClassId& id, std::istream& in);
 
    protected:
        inline virtual bool isEqual(const AbsDistributionND& other) const
        {
            const EllipticalKotz& r = static_cast<const EllipticalKotz&>(other);
            return EllipticalDistribution::isEqual(r) &&
                N_ == r.N_ && r_ == r.r_ && s_ == r.s_;
        }
 
    private:
        void calcNorm();
 
        inline EllipticalKotz(const EllipticalDistribution& b,
                              const unsigned N, const double r, const double s)
            : EllipticalDistribution(b),
              N_(N), r_(r), s_(s) {calcNorm();}
 
        unsigned N_;
        double r_;
        double s_;
        double norm_;
 
#ifdef SWIG
    public:
        inline EllipticalKotz(const double* location1, unsigned dim1,
                              const double* transform, unsigned nrows, unsigned ncols,
                              const unsigned N, const double r, const double s)
            : EllipticalDistribution(
                location1, dim1, transform, nrows, ncols,
                TransformedDistribution1D(PowerTransform1D(r, s), Gamma1D(0.0, 1.0, N/s)),
                TransformedDistribution1D(PowerTransform1D(r, s), Gamma1D(0.0, 1.0, (N+dim1/2.0-1)/s))),
                N_(N), r_(r), s_(s) {calcNorm();}
#endif // SWIG
    };
 
    class EllipticalPearsonTypeVII : public EllipticalDistribution
    {
    public:
        /**
        // In the constructor we must have s > 0, 2*N > dim1.
        // Also, the current implementation is limited to N < 7.
        */
        inline EllipticalPearsonTypeVII(const double* location1, unsigned dim1,
                                        const Matrix<double>& transformationMatrix,
                                        const unsigned N, const double s)
            : EllipticalDistribution(
                location1, dim1, transformationMatrix,
                PowerRatio1D(0.0, s, 0, N),
                PowerRatio1D(0.0, s, static_cast<int>(dim1)-2, N)) {}
 
        inline virtual ~EllipticalPearsonTypeVII() {}
 
        // Methods needed for I/O
        inline virtual gs::ClassId classId() const {return gs::ClassId(*this);}
        inline virtual bool write(std::ostream& os) const
        {
            return EllipticalDistribution::classId().write(os) &&
                   EllipticalDistribution::write(os);
        }
 
        static inline const char* classname()
            {return "npstat::EllipticalPearsonTypeVII";}
        static inline unsigned version() {return 1;}
        static inline EllipticalPearsonTypeVII* read(
            const gs::ClassId& id, std::istream& in)
        {
            static const gs::ClassId current(
                gs::ClassId::makeId<EllipticalPearsonTypeVII>());
            current.ensureSameId(id);
            CPP11_auto_ptr<EllipticalDistribution> pb =
                gs::read_obj<EllipticalDistribution>(in);
            return new EllipticalPearsonTypeVII(*pb);
        }
 
    private:
        inline EllipticalPearsonTypeVII(const EllipticalDistribution& b)
            : EllipticalDistribution(b) {}
 
#ifdef SWIG
    public:
        inline EllipticalPearsonTypeVII(const double* location1, unsigned dim1,
                                        const double* transform, unsigned nrows, unsigned ncols,
                                        const unsigned N, const double s)
            : EllipticalDistribution(
                location1, dim1, transform, nrows, ncols,
                PowerRatio1D(0.0, s, 0, N),
                PowerRatio1D(0.0, s, static_cast<int>(dim1)-2, N)) {}
#endif // SWIG
    };
 
    class EllipticalPearsonTypeII : public EllipticalDistribution
    {
    public:
        /**
        // Not sure why Pearson type II elliptical distribution does not have
        // the "s" parameter in the Fang's book. Here it is included, similar
        // to Pearson type VII.
        */
        inline EllipticalPearsonTypeII(const double* location1, unsigned dim1,
                                       const Matrix<double>& transformationMatrix,
                                       const unsigned m, const double s)
            : EllipticalDistribution(location1, dim1,
                                     transformationMatrix,
                                     Beta1D(0.0, s, 1.0, m+1U),
                                     Beta1D(0.0, s, dim1/2.0, m+1U)) {}
 
        inline virtual ~EllipticalPearsonTypeII() {}
 
        // Methods needed for I/O
        inline virtual gs::ClassId classId() const {return gs::ClassId(*this);}
        inline virtual bool write(std::ostream& os) const
        {
            return EllipticalDistribution::classId().write(os) &&
                   EllipticalDistribution::write(os);
        }
 
        static inline const char* classname()
            {return "npstat::EllipticalPearsonTypeII";}
        static inline unsigned version() {return 1;}
        static inline EllipticalPearsonTypeII* read(
            const gs::ClassId& id, std::istream& in)
        {
            static const gs::ClassId current(
                gs::ClassId::makeId<EllipticalPearsonTypeII>());
            current.ensureSameId(id);
            CPP11_auto_ptr<EllipticalDistribution> pb =
                gs::read_obj<EllipticalDistribution>(in);
            return new EllipticalPearsonTypeII(*pb);
        }
 
    private:
        inline EllipticalPearsonTypeII(const EllipticalDistribution& b)
            : EllipticalDistribution(b) {}
 
#ifdef SWIG
    public:
        inline EllipticalPearsonTypeII(const double* location1, unsigned dim1,
                                       const double* transform, unsigned nrows, unsigned ncols,
                                       const unsigned m, const double s)
            : EllipticalDistribution(location1, dim1,
                                     transform, nrows, ncols,
                                     Beta1D(0.0, s, 1.0, m+1U),
                                     Beta1D(0.0, s, dim1/2.0, m+1U)) {}
#endif // SWIG
    };
}
 
#endif // NPSTAT_ELLIPTICALDISTRIBUTIONS_HH_