exp_family.h

Go to the documentation of this file.

#ifndef EF_H
#define EF_H


#include "../shared_ptr.h"
#include "../base/bdmbase.h"
#include "../math/chmat.h"

namespace bdm
{


extern Uniform_RNG UniRNG;
extern Normal_RNG NorRNG;
extern Gamma_RNG GamRNG;

class eEF : public epdf
{
        public:
//      eEF() :epdf() {};
                eEF () : epdf () {};
                virtual double lognc() const = 0;

                virtual double evallog_nn (const vec &val) const {
                        bdm_error ("Not implemented");
                        return 0.0;
                }

                virtual double evallog (const vec &val) const {
                        double tmp;
                        tmp = evallog_nn (val) - lognc();
                        return tmp;
                }
                virtual vec evallog_m (const mat &Val) const {
                        vec x (Val.cols());
                        for (int i = 0;i < Val.cols();i++) {x (i) = evallog_nn (Val.get_col (i)) ;}
                        return x -lognc();
                }
                virtual vec evallog_m (const Array<vec> &Val) const {
                        vec x (Val.length());
                        for (int i = 0;i < Val.length();i++) {x (i) = evallog_nn (Val (i)) ;}
                        return x -lognc();
                }

                virtual void pow (double p) {
                        bdm_error ("Not implemented");
                }
};


class BMEF : public BM
{
        protected:
                double frg;
                double last_lognc;
        public:
                BMEF (double frg0 = 1.0) : BM (), frg (frg0) {}
                BMEF (const BMEF &B) : BM (B), frg (B.frg), last_lognc (B.last_lognc) {}
                virtual void set_statistics (const BMEF* BM0) {
                        bdm_error ("Not implemented");
                }

                virtual void bayes (const vec &data, const double w) {};
                //original Bayes
                void bayes (const vec &dt);

                virtual void flatten (const BMEF * B) {
                        bdm_error ("Not implemented");
                }

                BMEF* _copy_ () const {
                        bdm_error ("function _copy_ not implemented for this BM");
                        return NULL;
                }
};

template<class sq_T, template <typename> class TEpdf>
class mlnorm;

template<class sq_T>
class enorm : public eEF
{
        protected:
                vec mu;
                sq_T R;
        public:

                enorm () : eEF (), mu (), R () {};
                enorm (const vec &mu, const sq_T &R) {set_parameters (mu, R);}
                void set_parameters (const vec &mu, const sq_T &R);
                void from_setting (const Setting &root);
                void validate() {
                        bdm_assert_debug (mu.length() == R.rows(), "parameters mismatch");
                        dim = mu.length();
                }


                void dupdate (mat &v, double nu = 1.0);

                vec sample() const;

                double evallog_nn (const vec &val) const;
                double lognc () const;
                vec mean() const {return mu;}
                vec variance() const {return diag (R.to_mat());}
//      mlnorm<sq_T>* condition ( const RV &rvn ) const ; <=========== fails to cmpile. Why?
                shared_ptr<mpdf> condition ( const RV &rvn ) const;

                // target not typed to mlnorm<sq_T, enorm<sq_T> > &
                // because that doesn't compile (perhaps because we
                // haven't finished defining enorm yet), but the type
                // is required
                void condition ( const RV &rvn, mpdf &target ) const;

                shared_ptr<epdf> marginal (const RV &rvn ) const;
                void marginal ( const RV &rvn, enorm<sq_T> &target ) const;


                vec& _mu() {return mu;}
                const vec& _mu() const {return mu;}
                void set_mu (const vec mu0) { mu = mu0;}
                sq_T& _R() {return R;}
                const sq_T& _R() const {return R;}

};
UIREGISTER2 (enorm, chmat);
SHAREDPTR2 ( enorm, chmat );
UIREGISTER2 (enorm, ldmat);
SHAREDPTR2 ( enorm, ldmat );
UIREGISTER2 (enorm, fsqmat);
SHAREDPTR2 ( enorm, fsqmat );


class egiw : public eEF
{
        protected:
                ldmat V;
                double nu;
                int dimx;
                int nPsi;
        public:
                egiw() : eEF() {};
                egiw (int dimx0, ldmat V0, double nu0 = -1.0) : eEF() {set_parameters (dimx0, V0, nu0);};

                void set_parameters (int dimx0, ldmat V0, double nu0 = -1.0) {
                        dimx = dimx0;
                        nPsi = V0.rows() - dimx;
                        dim = dimx * (dimx + nPsi); // size(R) + size(Theta)

                        V = V0;
                        if (nu0 < 0) {
                                nu = 0.1 + nPsi + 2 * dimx + 2; // +2 assures finite expected value of R
                                // terms before that are sufficient for finite normalization
                        } else {
                                nu = nu0;
                        }
                }

                vec sample() const;
                vec mean() const;
                vec variance() const;

                vec est_theta() const;

                ldmat est_theta_cov() const;

                void mean_mat (mat &M, mat&R) const;
                double evallog_nn (const vec &val) const;
                double lognc () const;
                void pow (double p) {V *= p;nu *= p;};


                ldmat& _V() {return V;}
                const ldmat& _V() const {return V;}
                double& _nu()  {return nu;}
                const double& _nu() const {return nu;}
                void from_setting (const Setting &set) {
                        UI::get (nu, set, "nu", UI::compulsory);
                        UI::get (dimx, set, "dimx", UI::compulsory);
                        mat V;
                        UI::get (V, set, "V", UI::compulsory);
                        set_parameters (dimx, V, nu);
                        shared_ptr<RV> rv = UI::build<RV> (set, "rv", UI::compulsory);
                        set_rv (*rv);
                }
};
UIREGISTER ( egiw );
SHAREDPTR ( egiw );

class eDirich: public eEF
{
        protected:
                vec beta;
        public:

                eDirich () : eEF () {};
                eDirich (const eDirich &D0) : eEF () {set_parameters (D0.beta);};
                eDirich (const vec &beta0) {set_parameters (beta0);};
                void set_parameters (const vec &beta0) {
                        beta = beta0;
                        dim = beta.length();
                }

                vec sample() const {
                        bdm_error ("Not implemented");
                        return vec();
                }

                vec mean() const {return beta / sum (beta);};
                vec variance() const {double gamma = sum (beta); return elem_mult (beta, (beta + 1)) / (gamma* (gamma + 1));}
                double evallog_nn (const vec &val) const {
                        double tmp; tmp = (beta - 1) * log (val);
                        return tmp;
                }

                double lognc () const {
                        double tmp;
                        double gam = sum (beta);
                        double lgb = 0.0;
                        for (int i = 0;i < beta.length();i++) {lgb += lgamma (beta (i));}
                        tmp = lgb - lgamma (gam);
                        return tmp;
                }

                vec& _beta()  {return beta;}
};

class multiBM : public BMEF
{
        protected:
                eDirich est;
                vec &beta;
        public:
                multiBM () : BMEF (), est (), beta (est._beta()) {
                        if (beta.length() > 0) {last_lognc = est.lognc();}
                        else{last_lognc = 0.0;}
                }
                multiBM (const multiBM &B) : BMEF (B), est (B.est), beta (est._beta()) {}
                void set_statistics (const BM* mB0) {const multiBM* mB = dynamic_cast<const multiBM*> (mB0); beta = mB->beta;}
                void bayes (const vec &dt) {
                        if (frg < 1.0) {beta *= frg;last_lognc = est.lognc();}
                        beta += dt;
                        if (evalll) {ll = est.lognc() - last_lognc;}
                }
                double logpred (const vec &dt) const {
                        eDirich pred (est);
                        vec &beta = pred._beta();

                        double lll;
                        if (frg < 1.0)
                                {beta *= frg;lll = pred.lognc();}
                        else
                                if (evalll) {lll = last_lognc;}
                                else{lll = pred.lognc();}

                        beta += dt;
                        return pred.lognc() - lll;
                }
                void flatten (const BMEF* B) {
                        const multiBM* E = dynamic_cast<const multiBM*> (B);
                        // sum(beta) should be equal to sum(B.beta)
                        const vec &Eb = E->beta;//const_cast<multiBM*> ( E )->_beta();
                        beta *= (sum (Eb) / sum (beta));
                        if (evalll) {last_lognc = est.lognc();}
                }
                const eDirich& posterior() const {return est;};
                void set_parameters (const vec &beta0) {
                        est.set_parameters (beta0);
                        if (evalll) {last_lognc = est.lognc();}
                }
};

class egamma : public eEF
{
        protected:
                vec alpha;
                vec beta;
        public :
                egamma () : eEF (), alpha (0), beta (0) {};
                egamma (const vec &a, const vec &b) {set_parameters (a, b);};
                void set_parameters (const vec &a, const vec &b) {alpha = a, beta = b;dim = alpha.length();};

                vec sample() const;
                double evallog (const vec &val) const;
                double lognc () const;
                vec& _alpha() {return alpha;}
                vec& _beta() {return beta;}
                vec mean() const {return elem_div (alpha, beta);}
                vec variance() const {return elem_div (alpha, elem_mult (beta, beta)); }

                void from_setting (const Setting &set) {
                        epdf::from_setting (set); // reads rv
                        UI::get (alpha, set, "alpha", UI::compulsory);
                        UI::get (beta, set, "beta", UI::compulsory);
                        validate();
                }
                void validate() {
                        bdm_assert_debug (alpha.length() == beta.length(), "parameters do not match");
                        dim = alpha.length();
                }
};
UIREGISTER (egamma);
SHAREDPTR ( egamma );

class eigamma : public egamma
{
        protected:
        public :

                vec sample() const {return 1.0 / egamma::sample();};
                vec mean() const {return elem_div (beta, alpha - 1);}
                vec variance() const {vec mea = mean(); return elem_div (elem_mult (mea, mea), alpha - 2);}
};
/*
class emix : public epdf {
protected:
        int n;
        vec &w;
        Array<epdf*> Coms;
public:
        emix ( const RV &rv, vec &w0): epdf(rv), n(w0.length()), w(w0), Coms(n) {};
        void set_parameters( int &i, double wi, epdf* ep){w(i)=wi;Coms(i)=ep;}
        vec mean(){vec pom; for(int i=0;i<n;i++){pom+=Coms(i)->mean()*w(i);} return pom;};
};
*/


class euni: public epdf
{
        protected:
                vec low;
                vec high;
                vec distance;
                double nk;
                double lnk;
        public:
                euni () : epdf () {}
                euni (const vec &low0, const vec &high0) {set_parameters (low0, high0);}
                void set_parameters (const vec &low0, const vec &high0) {
                        distance = high0 - low0;
                        low = low0;
                        high = high0;
                        nk = prod (1.0 / distance);
                        lnk = log (nk);
                        dim = low.length();
                }

                double evallog (const vec &val) const  {
                        if (any (val < low) && any (val > high)) {return inf;}
                        else return lnk;
                }
                vec sample() const {
                        vec smp (dim);
#pragma omp critical
                        UniRNG.sample_vector (dim , smp);
                        return low + elem_mult (distance, smp);
                }
                vec mean() const {return (high -low) / 2.0;}
                vec variance() const {return (pow (high, 2) + pow (low, 2) + elem_mult (high, low)) / 3.0;}
                void from_setting (const Setting &set) {
                        epdf::from_setting (set); // reads rv and rvc

                        UI::get (high, set, "high", UI::compulsory);
                        UI::get (low, set, "low", UI::compulsory);
                        set_parameters(low,high);
                        validate();
                }
                void validate() {
                        bdm_assert(high.length()==low.length(), "Incompatible high and low vectors");
                        dim = high.length();
                        bdm_assert_debug (min (distance) > 0.0, "bad support");
                }
};
UIREGISTER(euni);

template < class sq_T, template <typename> class TEpdf = enorm >
class mlnorm : public mpdf_internal< TEpdf<sq_T> >
{
        protected:
                mat A;
                vec mu_const;
//                      vec& _mu; //cached epdf.mu; !!!!!! WHY NOT?
        public:
                mlnorm() : mpdf_internal< TEpdf<sq_T> >() {};
                mlnorm (const mat &A, const vec &mu0, const sq_T &R) : mpdf_internal< TEpdf<sq_T> >() {
                        set_parameters (A, mu0, R);
                }

                void set_parameters (const  mat &A0, const vec &mu0, const sq_T &R0) {
                        bdm_assert_debug (A0.rows() == mu0.length(), "mlnorm: A vs. mu mismatch");
                        bdm_assert_debug (A0.rows() == R0.rows(), "mlnorm: A vs. R mismatch");
                        
                        this->iepdf.set_parameters (zeros (A0.rows()), R0);
                        A = A0;
                        mu_const = mu0;
                        this->dimc = A0.cols();
                }
                void condition (const vec &cond) {
                        this->iepdf._mu() = A * cond + mu_const;
//R is already assigned;
                }

                const vec& _mu_const() const {return mu_const;}
                const mat& _A() const {return A;}
                mat _R() const { return this->iepdf._R().to_mat(); }

                template<typename sq_M>
                friend std::ostream &operator<< (std::ostream &os,  mlnorm<sq_M, enorm> &ml);

                void from_setting (const Setting &set) {
                        mpdf::from_setting (set);

                        UI::get (A, set, "A", UI::compulsory);
                        UI::get (mu_const, set, "const", UI::compulsory);
                        mat R0;
                        UI::get (R0, set, "R", UI::compulsory);
                        set_parameters (A, mu_const, R0);
                };
};
UIREGISTER2 (mlnorm,ldmat);
SHAREDPTR2 ( mlnorm, ldmat );
UIREGISTER2 (mlnorm,fsqmat);
SHAREDPTR2 ( mlnorm, fsqmat );
UIREGISTER2 (mlnorm, chmat);
SHAREDPTR2 ( mlnorm, chmat );

template<class sq_T>
class mgnorm : public mpdf_internal< enorm< sq_T > >
{
        private:
//                      vec &mu; WHY NOT?
                shared_ptr<fnc> g;

        public:
                mgnorm() : mpdf_internal<enorm<sq_T> >() { }
                inline void set_parameters (const shared_ptr<fnc> &g0, const sq_T &R0);
                inline void condition (const vec &cond);


                void from_setting (const Setting &set) {
                        shared_ptr<fnc> g = UI::build<fnc> (set, "g", UI::compulsory);

                        mat R;
                        vec dR;
                        if (UI::get (dR, set, "dR"))
                                R = diag (dR);
                        else
                                UI::get (R, set, "R", UI::compulsory);

                        set_parameters (g, R);
                }
};

UIREGISTER2 (mgnorm, chmat);
SHAREDPTR2 ( mgnorm, chmat );


class mlstudent : public mlnorm<ldmat, enorm>
{
        protected:
                ldmat Lambda;
                ldmat &_R;
                ldmat Re;
        public:
                mlstudent () : mlnorm<ldmat, enorm> (),
                                Lambda (),      _R (iepdf._R()) {}
                void set_parameters (const mat &A0, const vec &mu0, const ldmat &R0, const ldmat& Lambda0) {
                        iepdf.set_parameters (mu0, R0);// was Lambda, why?
                        A = A0;
                        mu_const = mu0;
                        Re = R0;
                        Lambda = Lambda0;
                }
                void condition (const vec &cond) {
                        iepdf._mu() = A * cond + mu_const;
                        double zeta;
                        //ugly hack!
                        if ( (cond.length() + 1) == Lambda.rows()) {
                                zeta = Lambda.invqform (concat (cond, vec_1 (1.0)));
                        } else {
                                zeta = Lambda.invqform (cond);
                        }
                        _R = Re;
                        _R *= (1 + zeta);// / ( nu ); << nu is in Re!!!!!!
                };

                void validate() {
                        bdm_assert_debug (A.rows() == mu_const.length(), "mlstudent: A vs. mu mismatch");
                        bdm_assert_debug (_R.rows() == A.rows(), "mlstudent: A vs. R mismatch");
                        
                }
};
class mgamma : public mpdf_internal<egamma>
{
        protected:

                double k;

                vec &_beta;

        public:
                mgamma() : mpdf_internal<egamma>(), k (0),
                                _beta (iepdf._beta()) {
                }

                void set_parameters (double k, const vec &beta0);

                void condition (const vec &val) {_beta = k / val;};
                void from_setting (const Setting &set) {
                        mpdf::from_setting (set); // reads rv and rvc
                        vec betatmp; // ugly but necessary
                        UI::get (betatmp, set, "beta", UI::compulsory);
                        UI::get (k, set, "k", UI::compulsory);
                        set_parameters (k, betatmp);
                }
};
UIREGISTER (mgamma);
SHAREDPTR (mgamma);

class migamma : public mpdf_internal<eigamma>
{
        protected:
                double k;

                vec &_alpha;

                vec &_beta;

        public:
                migamma() : mpdf_internal<eigamma>(),
                                k (0),
                                _alpha (iepdf._alpha()),
                                _beta (iepdf._beta()) {
                }

                migamma (const migamma &m) : mpdf_internal<eigamma>(),
                                k (0),
                                _alpha (iepdf._alpha()),
                                _beta (iepdf._beta()) {
                }

                void set_parameters (int len, double k0) {
                        k = k0;
                        iepdf.set_parameters ( (1.0 / (k*k) + 2.0) *ones (len) /*alpha*/, ones (len) /*beta*/);
                        dimc = dimension();
                };
                void condition (const vec &val) {
                        _beta = elem_mult (val, (_alpha - 1.0));
                };
};


class mgamma_fix : public mgamma
{
        protected:
                double l;
                vec refl;
        public:
                mgamma_fix () : mgamma (), refl () {};
                void set_parameters (double k0 , vec ref0, double l0) {
                        mgamma::set_parameters (k0, ref0);
                        refl = pow (ref0, 1.0 - l0);l = l0;
                        dimc = dimension();
                };

                void condition (const vec &val) {vec mean = elem_mult (refl, pow (val, l)); _beta = k / mean;};
};


class migamma_ref : public migamma
{
        protected:
                double l;
                vec refl;
        public:
                migamma_ref () : migamma (), refl () {};
                void set_parameters (double k0 , vec ref0, double l0) {
                        migamma::set_parameters (ref0.length(), k0);
                        refl = pow (ref0, 1.0 - l0);
                        l = l0;
                        dimc = dimension();
                };

                void condition (const vec &val) {
                        vec mean = elem_mult (refl, pow (val, l));
                        migamma::condition (mean);
                };

                void from_setting (const Setting &set);

                // TODO dodelat void to_setting( Setting &set ) const;
};


UIREGISTER (migamma_ref);
SHAREDPTR (migamma_ref);

class elognorm: public enorm<ldmat>
{
        public:
                vec sample() const {return exp (enorm<ldmat>::sample());};
                vec mean() const {vec var = enorm<ldmat>::variance();return exp (mu - 0.5*var);};

};

class mlognorm : public mpdf_internal<elognorm>
{
        protected:
                double sig2;

                vec &mu;
        public:
                mlognorm() : mpdf_internal<elognorm>(),
                                sig2 (0),
                                mu (iepdf._mu()) {
                }

                void set_parameters (int size, double k) {
                        sig2 = 0.5 * log (k * k + 1);
                        iepdf.set_parameters (zeros (size), 2*sig2*eye (size));

                        dimc = size;
                };

                void condition (const vec &val) {
                        mu = log (val) - sig2;//elem_mult ( refl,pow ( val,l ) );
                };

                void from_setting (const Setting &set);

                // TODO dodelat void to_setting( Setting &set ) const;

};

UIREGISTER (mlognorm);
SHAREDPTR (mlognorm);

class eWishartCh : public epdf
{
        protected:
                chmat Y;
                int p;
                double delta;
        public:
                void set_parameters (const mat &Y0, const double delta0) {Y = chmat (Y0);delta = delta0; p = Y.rows(); dim = p * p; }
                mat sample_mat() const {
                        mat X = zeros (p, p);

                        //sample diagonal
                        for (int i = 0;i < p;i++) {
                                GamRNG.setup (0.5* (delta - i) , 0.5);   // no +1 !! index if from 0
#pragma omp critical
                                X (i, i) = sqrt (GamRNG());
                        }
                        //do the rest
                        for (int i = 0;i < p;i++) {
                                for (int j = i + 1;j < p;j++) {
#pragma omp critical
                                        X (i, j) = NorRNG.sample();
                                }
                        }
                        return X*Y._Ch();// return upper triangular part of the decomposition
                }
                vec sample () const {
                        return vec (sample_mat()._data(), p*p);
                }
                void setY (const mat &Ch0) {copy_vector (dim, Ch0._data(), Y._Ch()._data());}
                void _setY (const vec &ch0) {copy_vector (dim, ch0._data(), Y._Ch()._data()); }
                const chmat& getY() const {return Y;}
};


class eiWishartCh: public epdf
{
        protected:
                eWishartCh W;
                int p;
                double delta;
        public:
                void set_parameters (const mat &Y0, const double delta0) {
                        delta = delta0;
                        W.set_parameters (inv (Y0), delta0);
                        dim = W.dimension(); p = Y0.rows();
                }
                vec sample() const {mat iCh; iCh = inv (W.sample_mat()); return vec (iCh._data(), dim);}
                void _setY (const vec &y0) {
                        mat Ch (p, p);
                        mat iCh (p, p);
                        copy_vector (dim, y0._data(), Ch._data());

                        iCh = inv (Ch);
                        W.setY (iCh);
                }
                virtual double evallog (const vec &val) const {
                        chmat X (p);
                        const chmat& Y = W.getY();

                        copy_vector (p*p, val._data(), X._Ch()._data());
                        chmat iX (p);X.inv (iX);
                        // compute
//                              \frac{ |\Psi|^{m/2}|X|^{-(m+p+1)/2}e^{-tr(\Psi X^{-1})/2} }{ 2^{mp/2}\Gamma_p(m/2)},
                        mat M = Y.to_mat() * iX.to_mat();

                        double log1 = 0.5 * p * (2 * Y.logdet()) - 0.5 * (delta + p + 1) * (2 * X.logdet()) - 0.5 * trace (M);
                        //Fixme! Multivariate gamma omitted!! it is ok for sampling, but not otherwise!!

                        /*                              if (0) {
                                                                mat XX=X.to_mat();
                                                                mat YY=Y.to_mat();

                                                                double log2 = 0.5*p*log(det(YY))-0.5*(delta+p+1)*log(det(XX))-0.5*trace(YY*inv(XX));
                                                                cout << log1 << "," << log2 << endl;
                                                        }*/
                        return log1;
                };

};

class rwiWishartCh : public mpdf_internal<eiWishartCh>
{
        protected:
                double sqd;
                vec refl;
                double l;
                int p;

        public:
                rwiWishartCh() : sqd (0), l (0), p (0) {}
                void set_parameters (int p0, double k, vec ref0, double l0) {
                        p = p0;
                        double delta = 2 / (k * k) + p + 3;
                        sqd = sqrt (delta - p - 1);
                        l = l0;
                        refl = pow (ref0, 1 - l);

                        iepdf.set_parameters (eye (p), delta);
                        dimc = iepdf.dimension();
                }
                void condition (const vec &c) {
                        vec z = c;
                        int ri = 0;
                        for (int i = 0;i < p*p;i += (p + 1)) {//trace diagonal element
                                z (i) = pow (z (i), l) * refl (ri);
                                ri++;
                        }

                        iepdf._setY (sqd*z);
                }
};

enum RESAMPLING_METHOD { MULTINOMIAL = 0, STRATIFIED = 1, SYSTEMATIC = 3 };
class eEmp: public epdf
{
        protected :
                int n;
                vec w;
                Array<vec> samples;
        public:
                eEmp () : epdf (), w (), samples () {};
                eEmp (const eEmp &e) : epdf (e), w (e.w), samples (e.samples) {};

                void set_statistics (const vec &w0, const epdf &pdf0);
                void set_statistics (const epdf &pdf0 , int n) {set_statistics (ones (n) / n, pdf0);};
                void set_samples (const epdf* pdf0);
                void set_parameters (int n0, bool copy = true) {n = n0; w.set_size (n0, copy);samples.set_size (n0, copy);};
                void set_parameters (const Array<vec> &Av) {
                        bdm_assert_debug(Av.size()>0,"Empty samples"); 
                        n = Av.size(); 
                        epdf::set_parameters(Av(0).length());
                        w=1/n*ones(n);
                        samples=Av;
                };
                vec& _w()  {return w;};
                const vec& _w() const {return w;};
                Array<vec>& _samples() {return samples;};
                const Array<vec>& _samples() const {return samples;};
                ivec resample (RESAMPLING_METHOD method = SYSTEMATIC);

                vec sample() const {
                        bdm_error ("Not implemented");
                        return vec();
                }

                double evallog (const vec &val) const {
                        bdm_error ("Not implemented");
                        return 0.0;
                }

                vec mean() const {
                        vec pom = zeros (dim);
                        for (int i = 0;i < n;i++) {pom += samples (i) * w (i);}
                        return pom;
                }
                vec variance() const {
                        vec pom = zeros (dim);
                        for (int i = 0;i < n;i++) {pom += pow (samples (i), 2) * w (i);}
                        return pom -pow (mean(), 2);
                }
                void qbounds (vec &lb, vec &ub, double perc = 0.95) const {
                        // lb in inf so than it will be pushed below;
                        lb.set_size (dim);
                        ub.set_size (dim);
                        lb = std::numeric_limits<double>::infinity();
                        ub = -std::numeric_limits<double>::infinity();
                        int j;
                        for (int i = 0;i < n;i++) {
                                for (j = 0;j < dim; j++) {
                                        if (samples (i) (j) < lb (j)) {lb (j) = samples (i) (j);}
                                        if (samples (i) (j) > ub (j)) {ub (j) = samples (i) (j);}
                                }
                        }
                }
};



template<class sq_T>
void enorm<sq_T>::set_parameters (const vec &mu0, const sq_T &R0)
{
//Fixme test dimensions of mu0 and R0;
        mu = mu0;
        R = R0;
        validate();
};

template<class sq_T>
void enorm<sq_T>::from_setting (const Setting &set)
{
        epdf::from_setting (set); //reads rv

        UI::get (mu, set, "mu", UI::compulsory);
        mat Rtmp;// necessary for conversion
        UI::get (Rtmp, set, "R", UI::compulsory);
        R = Rtmp; // conversion
        validate();
}

template<class sq_T>
void enorm<sq_T>::dupdate (mat &v, double nu)
{
        //
};

// template<class sq_T>
// void enorm<sq_T>::tupdate ( double phi, mat &vbar, double nubar ) {
//      //
// };

template<class sq_T>
vec enorm<sq_T>::sample() const
{
        vec x (dim);
#pragma omp critical
        NorRNG.sample_vector (dim, x);
        vec smp = R.sqrt_mult (x);

        smp += mu;
        return smp;
};

// template<class sq_T>
// double enorm<sq_T>::eval ( const vec &val ) const {
//      double pdfl,e;
//      pdfl = evallog ( val );
//      e = exp ( pdfl );
//      return e;
// };

template<class sq_T>
double enorm<sq_T>::evallog_nn (const vec &val) const
{
        // 1.83787706640935 = log(2pi)
        double tmp = -0.5 * (R.invqform (mu - val));// - lognc();
        return  tmp;
};

template<class sq_T>
inline double enorm<sq_T>::lognc () const
{
        // 1.83787706640935 = log(2pi)
        double tmp = 0.5 * (R.cols() * 1.83787706640935 + R.logdet());
        return tmp;
};


// template<class sq_T>
// vec mlnorm<sq_T>::samplecond (const  vec &cond, double &lik ) {
//      this->condition ( cond );
//      vec smp = epdf.sample();
//      lik = epdf.eval ( smp );
//      return smp;
// }

// template<class sq_T>
// mat mlnorm<sq_T>::samplecond (const vec &cond, vec &lik, int n ) {
//      int i;
//      int dim = rv.count();
//      mat Smp ( dim,n );
//      vec smp ( dim );
//      this->condition ( cond );
//
//      for ( i=0; i<n; i++ ) {
//              smp = epdf.sample();
//              lik ( i ) = epdf.eval ( smp );
//              Smp.set_col ( i ,smp );
//      }
//
//      return Smp;
// }


template<class sq_T>
shared_ptr<epdf> enorm<sq_T>::marginal ( const RV &rvn ) const
{
        enorm<sq_T> *tmp = new enorm<sq_T> ();
        shared_ptr<epdf> narrow(tmp);
        marginal ( rvn, *tmp );
        return narrow;
}

template<class sq_T>
void enorm<sq_T>::marginal ( const RV &rvn, enorm<sq_T> &target ) const
{
        bdm_assert_debug (isnamed(), "rv description is not assigned");
        ivec irvn = rvn.dataind (rv);

        sq_T Rn (R, irvn);  // select rows and columns of R

        target.set_rv ( rvn );
        target.set_parameters (mu (irvn), Rn);
}

template<class sq_T>
shared_ptr<mpdf> enorm<sq_T>::condition ( const RV &rvn ) const
{
        mlnorm<sq_T> *tmp = new mlnorm<sq_T> ();
        shared_ptr<mpdf> narrow(tmp);
        condition ( rvn, *tmp );
        return narrow;
}

template<class sq_T>
void enorm<sq_T>::condition ( const RV &rvn, mpdf &target ) const
{
        typedef mlnorm<sq_T> TMlnorm;

        bdm_assert_debug (isnamed(), "rvs are not assigned");
        TMlnorm &uptarget = dynamic_cast<TMlnorm &>(target);

        RV rvc = rv.subt (rvn);
        bdm_assert_debug ( (rvc._dsize() + rvn._dsize() == rv._dsize()), "wrong rvn");
        //Permutation vector of the new R
        ivec irvn = rvn.dataind (rv);
        ivec irvc = rvc.dataind (rv);
        ivec perm = concat (irvn , irvc);
        sq_T Rn (R, perm);

        //fixme - could this be done in general for all sq_T?
        mat S = Rn.to_mat();
        //fixme
        int n = rvn._dsize() - 1;
        int end = R.rows() - 1;
        mat S11 = S.get (0, n, 0, n);
        mat S12 = S.get (0, n , rvn._dsize(), end);
        mat S22 = S.get (rvn._dsize(), end, rvn._dsize(), end);

        vec mu1 = mu (irvn);
        vec mu2 = mu (irvc);
        mat A = S12 * inv (S22);
        sq_T R_n (S11 - A *S12.T());

        uptarget.set_rv (rvn);
        uptarget.set_rvc (rvc);
        uptarget.set_parameters (A, mu1 - A*mu2, R_n);
}

template<class sq_T>
void mgnorm<sq_T >::set_parameters (const shared_ptr<fnc> &g0, const sq_T &R0) {
        g = g0;
        this->iepdf.set_parameters (zeros (g->dimension()), R0);
}

template<class sq_T>
void mgnorm<sq_T >::condition (const vec &cond) {this->iepdf._mu() = g->eval (cond);};

template<class sq_T>
std::ostream &operator<< (std::ostream &os,  mlnorm<sq_T> &ml)
{
        os << "A:" << ml.A << endl;
        os << "mu:" << ml.mu_const << endl;
        os << "R:" << ml._R() << endl;
        return os;
};

}
#endif //EF_H

---