root/bdm/stat/emix.h @ 229

/*!
  \file
  \brief Probability distributions for Mixtures of pdfs
  \author Vaclav Smidl.

  -----------------------------------
  BDM++ - C++ library for Bayesian Decision Making under Uncertainty

  Using IT++ for numerical operations
  -----------------------------------
*/

#ifndef MX_H
#define MX_H

#include "libBM.h"
#include "libEF.h"
//#include <std>

using namespace itpp;

//this comes first because it is used inside emix!

/*! \brief Class representing ratio of two densities
which arise e.g. by applying the Bayes rule.
It represents density in the form:
\f[
f(rv|rvc) = \frac{f(rv,rvc)}{f(rvc)}
\f]
where \f$ f(rvc) = \int f(rv,rvc) d\ rv \f$.

In particular this type of arise by conditioning of a mixture model.

At present the only supported operation is evallogcond().
 */
class mratio: public mpdf {
protected:
        //! Nominator in the form of mpdf
        const epdf* nom;
        //!Denominator in the form of epdf
        epdf* den;
        //!flag for destructor
        bool destroynom;
        //!datalink between conditional and nom
        datalink_m2e dl;
public:
        //!Default constructor. By default, the given epdf is not copied!
        //! It is assumed that this function will be used only temporarily.
        mratio ( const epdf* nom0, const RV &rv, bool copy=false ) :mpdf ( rv,nom0->_rv().subt ( rv ) ), dl ( rv,rvc,nom0->_rv() ) {
                if ( copy ) {
//                      nom = nom0->_copy_();
                        it_error ( "todo" );
                        destroynom=true;
                }
                else {
                        nom = nom0;
                        destroynom = false;
                }
                it_assert_debug ( rvc.length() >0,"Makes no sense to use this object!" );
                den = nom->marginal ( rvc );
        };
        double evallogcond ( const vec &val, const vec &cond ) {
                double tmp;
                vec nom_val ( rv.count() +rvc.count() );
                dl.fill_val_cond ( nom_val,val,cond );
                tmp = exp ( nom->evallog ( nom_val ) - den->evallog ( cond ) );
                it_assert_debug(std::isfinite(tmp),"Infinite value");
                return tmp;
        }
        //! Object takes ownership of nom and will destroy it
        void ownnom() {destroynom=true;}
        //! Default destructor
        ~mratio() {delete den; if ( destroynom ) {delete nom;}}
};

/*!
* \brief Mixture of epdfs

Density function:
\f[
f(x) = \sum_{i=1}^{n} w_{i} f_i(x), \quad \sum_{i=1}^n w_i = 1.
\f]
where \f$f_i(x)\f$ is any density on random variable \f$x\f$, called \a component,

*/
class emix : public epdf {
protected:
        //! weights of the components
        vec w;
        //! Component (epdfs)
        Array<epdf*> Coms;
        //!Flag if owning Coms
        bool destroyComs;
public:
        //!Default constructor
        emix ( const RV &rv ) : epdf ( rv ) {};
        //! Set weights \c w and components \c Coms
        //!By default Coms are copied inside. Parameter \c copy can be set to false if Coms live externally. Use method ownComs() if Coms should be destroyed by the destructor.
        void set_parameters ( const vec &w, const Array<epdf*> &Coms, bool copy=true );

        vec sample() const;
        vec mean() const {
                int i; vec mu = zeros ( rv.count() );
                for ( i = 0;i < w.length();i++ ) {mu += w ( i ) * Coms ( i )->mean(); }
                return mu;
        }
        vec variance() const {
                //non-central moment
                vec mom2 = zeros(rv.count());
                for ( int i = 0;i < w.length();i++ ) {mom2 += w ( i ) * pow(Coms ( i )->mean(),2); }
                //central moment
                return mom2-pow(mean(),2);
        }
        double evallog ( const vec &val ) const {
                int i;
                double sum = 0.0;
                for ( i = 0;i < w.length();i++ ) {sum += w ( i ) * exp ( Coms ( i )->evallog ( val ) );}
                if (sum==0.0){sum=std::numeric_limits<double>::epsilon();}
                double tmp=log ( sum );
                it_assert_debug(std::isfinite(tmp),"Infinite");
                return tmp;
        };
        vec evallog_m ( const mat &Val ) const {
                vec x=zeros ( Val.cols() );
                for ( int i = 0; i < w.length(); i++ ) {
                        x+= w ( i ) *exp ( Coms ( i )->evallog_m ( Val ) );
                }
                return log ( x );
        };
        //! Auxiliary function that returns pdflog for each component
        mat evallog_M ( const mat &Val ) const {
                mat X ( w.length(), Val.cols() );
                for ( int i = 0; i < w.length(); i++ ) {
                        X.set_row ( i, w ( i ) *exp ( Coms ( i )->evallog_m ( Val ) ) );
                }
                return X;
        };

        emix* marginal ( const RV &rv ) const;
        mratio* condition ( const RV &rv ) const; //why not mratio!!

//Access methods
        //! returns a pointer to the internal mean value. Use with Care!
        vec& _w() {return w;}
        virtual ~emix() {if ( destroyComs ) {for ( int i=0;i<Coms.length();i++ ) {delete Coms ( i );}}}
        //! Auxiliary function for taking ownership of the Coms()
        void ownComs() {destroyComs=true;}

        //!access function
        epdf* _Coms ( int i ) {return Coms ( i );}
};

/*! \brief Chain rule decomposition of epdf

Probability density in the form of Chain-rule decomposition:
\[
f(x_1,x_2,x_3) = f(x_1|x_2,x_3)f(x_2,x_3)f(x_3)
\]
Note that
*/
class mprod: public compositepdf, public mpdf {
protected:
        //! pointers to epdfs - shortcut to mpdfs()._epdf()
        Array<epdf*> epdfs;
        //! Data link for each mpdfs
        Array<datalink_m2m*> dls;
public:
        /*!\brief Constructor from list of mFacs,
        */
        mprod ( Array<mpdf*> mFacs ) : compositepdf ( mFacs ), mpdf ( getrv ( true ),RV() ), epdfs ( n ), dls ( n ) {
                setrvc ( rv,rvc );
                // rv and rvc established = > we can link them with mpdfs
                for ( int i = 0;i < n;i++ ) {
                        dls ( i ) = new datalink_m2m ( mpdfs ( i )->_rv(), mpdfs ( i )->_rvc(), rv, rvc );
                }

                for ( int i=0;i<n;i++ ) {
                        epdfs ( i ) =& ( mpdfs ( i )->_epdf() );
                }
        };

        double evallogcond ( const vec &val, const vec &cond ) {
                int i;
                double res = 1.0;
                for ( i = n - 1;i >= 0;i-- ) {
                        /*                      if ( mpdfs(i)->_rvc().count() >0) {
                                                        mpdfs ( i )->condition ( dls ( i )->get_cond ( val,cond ) );
                                                }
                                                // add logarithms
                                                res += epdfs ( i )->evallog ( dls ( i )->get_val ( val ) );*/
                        res *= mpdfs ( i )->evallogcond (
                                   dls ( i )->get_val ( val ),
                                   dls ( i )->get_cond ( val, cond )
                               );
                }
                return res;
        }
        vec samplecond ( const vec &cond, double &ll ) {
                //! Ugly hack to help to discover if mpfs are not in proper order. Correct solution = check that explicitely.
                vec smp= std::numeric_limits<double>::infinity() * ones ( rv.count() );
                vec smpi;
                ll = 0;
                // Hard assumption here!!! We are going backwards, to assure that samples that are needed from smp are already generated!
                for ( int i = ( n - 1 );i >= 0;i-- ) {
                        if ( mpdfs ( i )->_rvc().count() ) {
                                mpdfs ( i )->condition ( dls ( i )->get_cond ( smp ,cond ) ); // smp is val here!!
                        }
                        smpi = epdfs ( i )->sample();
                        // copy contribution of this pdf into smp
                        dls ( i )->fill_val ( smp, smpi );
                        // add ith likelihood contribution
                        ll+=epdfs ( i )->evallog ( smpi );
                }
                return smp;
        }
        mat samplecond ( const vec &cond, vec &ll, int N ) {
                mat Smp ( rv.count(),N );
                for ( int i=0;i<N;i++ ) {Smp.set_col ( i,samplecond ( cond,ll ( i ) ) );}
                return Smp;
        }

        ~mprod() {};
};

//! Product of independent epdfs. For dependent pdfs, use mprod.
class eprod: public epdf {
protected:
        //! Components (epdfs)
        Array<const epdf*> epdfs;
        //! Array of indeces
        Array<datalink_e2e*> dls;
public:
        eprod ( const Array<const epdf*> epdfs0 ) : epdf ( RV() ),epdfs ( epdfs0 ),dls ( epdfs.length() ) {
                bool independent=true;
                for ( int i=0;i<epdfs.length();i++ ) {
                        independent=rv.add ( epdfs ( i )->_rv() );
                        it_assert_debug ( independent==true, "eprod:: given components are not independent ." );
                }
                for ( int i=0;i<epdfs.length();i++ ) {
                        dls ( i ) = new datalink_e2e ( epdfs ( i )->_rv() , rv );
                }
        }

        vec mean() const {
                vec tmp ( rv.count() );
                for ( int i=0;i<epdfs.length();i++ ) {
                        vec pom = epdfs ( i )->mean();
                        dls ( i )->fill_val ( tmp, pom );
                }
                return tmp;
        }
        vec variance() const {
                vec tmp ( rv.count() ); //second moment
                for ( int i=0;i<epdfs.length();i++ ) {
                        vec pom = epdfs ( i )->mean();
                        dls ( i )->fill_val ( tmp, pow(pom,2) );
                }
                return tmp-pow(mean(),2);
        }
        vec sample() const {
                vec tmp ( rv.count() );
                for ( int i=0;i<epdfs.length();i++ ) {
                        vec pom = epdfs ( i )->sample();
                        dls ( i )->fill_val ( tmp, pom );
                }
                return tmp;
        }
        double evallog ( const vec &val ) const {
                double tmp=0;
                for ( int i=0;i<epdfs.length();i++ ) {
                        tmp+=epdfs ( i )->evallog ( dls ( i )->get_val ( val ) );
                }
                it_assert_debug(std::isfinite(tmp),"Infinite");
                return tmp;
        }
        //!access function
        const epdf* operator () ( int i ) const {it_assert_debug ( i<epdfs.length(),"wrong index" );return epdfs ( i );}

        //!Destructor
        ~eprod() {for ( int i=0;i<epdfs.length();i++ ) {delete dls ( i );}}
};


/*! \brief Mixture of mpdfs with constant weights, all mpdfs are of equal type

*/
class mmix : public mpdf {
protected:
        //! Component (epdfs)
        Array<mpdf*> Coms;
        //!Internal epdf
        emix Epdf;
public:
        //!Default constructor
        mmix ( RV &rv, RV &rvc ) : mpdf ( rv, rvc ), Epdf ( rv ) {ep = &Epdf;};
        //! Set weights \c w and components \c R
        void set_parameters ( const vec &w, const Array<mpdf*> &Coms ) {
                Array<epdf*> Eps ( Coms.length() );

                for ( int i = 0;i < Coms.length();i++ ) {
                        Eps ( i ) = & ( Coms ( i )->_epdf() );
                }
                Epdf.set_parameters ( w, Eps );
        };

        void condition ( const vec &cond ) {
                for ( int i = 0;i < Coms.length();i++ ) {Coms ( i )->condition ( cond );}
        };
};

#endif //MX_H