root/library/bdm/stat/emix.h @ 1068

Revision 1068, 14.2 kB (checked in by mido, 14 years ago)

patch of documentation - all conditional pdfs revised

  • Property svn:eol-style set to native
RevLine 
[107]1/*!
2  \file
3  \brief Probability distributions for Mixtures of pdfs
4  \author Vaclav Smidl.
5
6  -----------------------------------
7  BDM++ - C++ library for Bayesian Decision Making under Uncertainty
8
9  Using IT++ for numerical operations
10  -----------------------------------
11*/
12
[394]13#ifndef EMIX_H
14#define EMIX_H
[107]15
[477]16#define LOG2  0.69314718055995
[333]17
[461]18#include "../shared_ptr.h"
[384]19#include "exp_family.h"
[107]20
[286]21namespace bdm {
[107]22
[182]23//this comes first because it is used inside emix!
24
25/*! \brief Class representing ratio of two densities
[1068]26
[182]27which arise e.g. by applying the Bayes rule.
28It represents density in the form:
29\f[
30f(rv|rvc) = \frac{f(rv,rvc)}{f(rvc)}
31\f]
32where \f$ f(rvc) = \int f(rv,rvc) d\ rv \f$.
33
34In particular this type of arise by conditioning of a mixture model.
35
[211]36At present the only supported operation is evallogcond().
[182]37 */
[693]38class mratio: public pdf {
[192]39protected:
[1064]40    //! Nominator in the form of pdf
41    const epdf* nom;
[504]42
[1064]43    //!Denominator in the form of epdf
44    shared_ptr<epdf> den;
[504]45
[1064]46    //!flag for destructor
47    bool destroynom;
48    //!datalink between conditional and nom
49    datalink_m2e dl;
[192]50public:
[1064]51    //!Default constructor. By default, the given epdf is not copied!
52    //! It is assumed that this function will be used only temporarily.
53    mratio ( const epdf* nom0, const RV &rv, bool copy = false ) : pdf ( ), dl ( ) {
54        // adjust rv and rvc
[675]55
[1064]56        set_rv ( rv );
57        dim = rv._dsize();
[675]58
[1064]59        rvc = nom0->_rv().subt ( rv );
60        dimc = rvc._dsize();
[737]61
[1064]62        //prepare data structures
63        if ( copy ) {
64            bdm_error ( "todo" );
65            // destroynom = true;
66        } else {
67            nom = nom0;
68            destroynom = false;
69        }
70        bdm_assert_debug ( rvc.length() > 0, "Makes no sense to use this object!" );
[477]71
[1064]72        // build denominator
73        den = nom->marginal ( rvc );
74        dl.set_connection ( rv, rvc, nom0->_rv() );
75    };
[675]76
[1064]77    double evallogcond ( const vec &val, const vec &cond ) {
78        double tmp;
79        vec nom_val ( dimension() + dimc );
80        dl.pushup_cond ( nom_val, val, cond );
81        tmp = exp ( nom->evallog ( nom_val ) - den->evallog ( cond ) );
82        return tmp;
83    }
[766]84
[1064]85    //! Returns a sample from the density conditioned on \c cond, \f$x \sim epdf(rv|cond)\f$. \param cond is numeric value of \c rv
86    virtual vec samplecond ( const vec &cond ) NOT_IMPLEMENTED(0);
[766]87
[1064]88    //! Object takes ownership of nom and will destroy it
89    void ownnom() {
90        destroynom = true;
91    }
92    //! Default destructor
93    ~mratio() {
94        if ( destroynom ) {
95            delete nom;
96        }
97    }
[550]98
[766]99
[550]100private:
[1064]101    // not implemented
102    mratio ( const mratio & );
103    mratio &operator= ( const mratio & );
[182]104};
105
[886]106class emix; //forward
107
[1066]108//! \brief Base class (interface) for mixtures
[886]109class emix_base : public epdf {
[1064]110protected:
111    //! reference to vector of weights
112    vec &w;
113    //! function returning ith component
114    virtual const epdf * component(const int &i) const=0;
115
116    virtual int no_coms() const=0;
117
118public:
119
120    emix_base(vec &w0): w(w0) {}
121
122    void validate ();
123
124    vec sample() const;
125
126    vec mean() const;
127
128    vec variance() const;
129
130    double evallog ( const vec &val ) const;
131
132    vec evallog_mat ( const mat &Val ) const;
133
134    //! Auxiliary function that returns pdflog for each component
135    mat evallog_coms ( const mat &Val ) const;
136
137    shared_ptr<epdf> marginal ( const RV &rv ) const;
138    //! Update already existing marginal density  \c target
139    void marginal ( const RV &rv, emix &target ) const;
140    shared_ptr<pdf> condition ( const RV &rv ) const;
141
142    //Access methods
143    //! returns a reference to the internal weights. Use with Care!
144    vec& _w() {
145        return w;
146    }
147
148    const vec& _w() const {
149        return w;
150    }
151    //!access
152    const epdf* _com(int i) const {
153        return component(i);
154    }
155
[886]156};
157
[107]158/*!
159* \brief Mixture of epdfs
160
161Density function:
162\f[
163f(x) = \sum_{i=1}^{n} w_{i} f_i(x), \quad \sum_{i=1}^n w_i = 1.
164\f]
165where \f$f_i(x)\f$ is any density on random variable \f$x\f$, called \a component,
166
167*/
[886]168class emix : public emix_base {
[162]169protected:
[1064]170    //! weights of the components
171    vec weights;
[559]172
[1064]173    //! Component (epdfs)
174    Array<shared_ptr<epdf> > Coms;
[504]175
[162]176public:
[1064]177    //! Default constructor
178    emix ( ) : emix_base ( weights) { }
[559]179
[1064]180    const epdf* component(const int &i) const {
181        return Coms(i).get();
182    }
183    void validate();
[107]184
185
[1064]186    int no_coms() const {
187        return Coms.length();
188    }
189
[1066]190    /*! Create object from the following structure
191
192    \code
193    class = 'emix';
194
195    pdfs = { list of any bdm::pdf offsprings };   % pdfs in the mixture, bdm::pdf::from_setting
196    weights = [... ];                             % vector of weights of pdfs in the mixture
197    --- inherited fields ---
198    bdm::emix_base::from_setting
199    \endcode
200
201    */
[1064]202    void from_setting ( const Setting &set );
[1066]203
[1064]204    void to_setting  (Setting  &set) const;
205
206    void set_rv ( const RV &rv ) {
207        epdf::set_rv ( rv );
208        for ( int i = 0; i < no_coms(); i++ ) {
209            Coms( i )->set_rv ( rv );
210        }
211    }
212
213    Array<shared_ptr<epdf> >& _Coms ( ) {
214        return Coms;
215    }
[333]216};
[886]217SHAREDPTR ( emix );
218UIREGISTER ( emix );
[333]219
[886]220
[115]221/*! \brief Chain rule decomposition of epdf
222
[145]223Probability density in the form of Chain-rule decomposition:
224\[
225f(x_1,x_2,x_3) = f(x_1|x_2,x_3)f(x_2,x_3)f(x_3)
226\]
227Note that
[115]228*/
[693]229class mprod: public pdf {
[507]230private:
[1064]231    Array<shared_ptr<pdf> > pdfs;
[507]232
[1064]233    //! Data link for each pdfs
234    Array<shared_ptr<datalink_m2m> > dls;
[461]235
[162]236public:
[1064]237    //! \brief Default constructor
238    mprod() { }
[507]239
[1064]240    /*!\brief Constructor from list of mFacs
241    */
242    mprod ( const Array<shared_ptr<pdf> > &mFacs ) {
243        set_elements ( mFacs );
244    }
245    //! Set internal \c pdfs from given values
246    void set_elements ( const Array<shared_ptr<pdf> > &mFacs );
[477]247
[1064]248    double evallogcond ( const vec &val, const vec &cond );
[395]249
[1064]250    vec evallogcond_mat ( const mat &Dt, const vec &cond );
[395]251
[1064]252    vec evallogcond_mat ( const Array<vec> &Dt, const vec &cond );
[739]253
[1064]254    //TODO smarter...
255    vec samplecond ( const vec &cond ) {
256        //! Ugly hack to help to discover if mpfs are not in proper order. Correct solution = check that explicitely.
257        vec smp = std::numeric_limits<double>::infinity() * ones ( dimension() );
258        vec smpi;
259        // Hard assumption here!!! We are going backwards, to assure that samples that are needed from smp are already generated!
260        for ( int i = ( pdfs.length() - 1 ); i >= 0; i-- ) {
261            // generate contribution of this pdf
262            smpi = pdfs ( i )->samplecond ( dls ( i )->get_cond ( smp , cond ) );
263            // copy contribution of this pdf into smp
264            dls ( i )->pushup ( smp, smpi );
265        }
266        return smp;
267    }
[162]268
[1068]269    //! Create object from the following structure
270    //! \code
271    //! class='mprod';
272    //! pdfs = { list of bdm::pdf };          % list of bdm::pdf offsprings stored in the order of chain rule, bdm::pdf::from_setting
273    //! --- inherited fields ---
274    //! bdm::pdf::from_setting
[1064]275    //! \endcode
276    //!@}
277    void from_setting ( const Setting &set ) ;
[1068]278    void to_setting  (Setting  &set) const;   
[107]279};
[477]280UIREGISTER ( mprod );
[529]281SHAREDPTR ( mprod );
[107]282
[886]283
[1066]284//! \brief Base class (interface) for bdm::eprod
[886]285class eprod_base: public epdf {
[168]286protected:
[1064]287    //! Array of indices
288    Array<datalink*> dls;
289    //! interface for a factor
[979]290public:
[1064]291    virtual const epdf* factor(int i) const NOT_IMPLEMENTED(NULL);
292    //!number of factors
293    virtual const int no_factors() const NOT_IMPLEMENTED(0);
294    //! Default constructor
295    eprod_base () :  dls (0) {};
296    //! Set internal
297    vec mean() const;
[886]298
[1064]299    vec variance() const;
[886]300
[1064]301    vec sample() const;
[886]302
[1064]303    double evallog ( const vec &val ) const;
[886]304
[1064]305    //!Destructor
306    ~eprod_base() {
307        for ( int i = 0; i < dls.length(); i++ ) {
308            delete dls ( i );
309        }
310    }
311    void validate() {
312        epdf::validate();
313        dls.set_length ( no_factors() );
314
315        bool independent = true;
316        dim = 0;
317        for ( int i = 0; i < no_factors(); i++ ) {
318            independent = rv.add ( factor ( i )->_rv() );
319            dim += factor ( i )->dimension();
320            bdm_assert_debug ( independent, "eprod:: given components are not independent." );
321        };
322
323        //
324        int cumdim = 0;
325        int dimi = 0;
326        int i;
327        for ( i = 0; i < no_factors(); i++ ) {
328            dls ( i ) = new datalink;
329            if ( isnamed() ) { // rvs are complete
330                dls ( i )->set_connection ( factor ( i )->_rv() , rv );
331            } else { //rvs are not reliable
332                dimi = factor ( i )->dimension();
333                dls ( i )->set_connection ( dimi, dim, linspace ( cumdim, cumdim + dimi - 1 ) );
334                cumdim += dimi;
335            }
336        }
337
338    }
[886]339};
[168]340
[1066]341//! \brief Product of independent epdfs. For dependent pdfs, use bdm::mprod.
[1064]342class eprod: public eprod_base {
343protected:
344    Array<shared_ptr<epdf> > factors;
345public:
346    const epdf* factor(int i) const {
347        return factors(i).get();
348    }
349    const int no_factors() const {
350        return factors.length();
351    }
352    void set_parameters ( const Array<shared_ptr<epdf> > &epdfs0) {
353        factors = epdfs0;
354    }
[1066]355
356    /*! Create object from the following structure
357
358    \code
359    class = 'eprod';
360    pdfs = { list of any bdm::epdf offsprings };   % pdfs in the product, bdm::epdf::from_setting
361    --- inherited fields ---
362    bdm::eprod_base::from_setting
363    \endcode
364
365    */
[1064]366    void from_setting(const Setting &set) {
367        UI::get(factors,set,"pdfs",UI::compulsory);
368    }
[886]369};
[944]370UIREGISTER(eprod);
[204]371
[1068]372//! \brief Internal class similar to eprod - factors are external pointers. To be used internally! 
[1064]373class eprod_internal: public eprod_base {
374protected:
375    Array<epdf* > factors;
376    const epdf* factor(int i) const {
377        return factors(i);
378    }
379    const int no_factors() const {
380        return factors.length();
381    }
382public:
383    void set_parameters ( const Array<epdf *> &epdfs0) {
384        factors = epdfs0;
385    }
[168]386};
387
[693]388/*! \brief Mixture of pdfs with constant weights, all pdfs are of equal RV and RVC
[124]389
390*/
[693]391class mmix : public pdf {
[488]392protected:
[1064]393    //! Component (pdfs)
394    Array<shared_ptr<pdf> > Coms;
395    //!weights of the components
396    vec w;
[488]397public:
[1064]398    //!Default constructor
399    mmix() : Coms ( 0 ) { };
[461]400
[1064]401    double evallogcond ( const vec &dt, const vec &cond ) {
402        double ll = 0.0;
403        for ( int i = 0; i < Coms.length(); i++ ) {
404            ll += Coms ( i )->evallogcond ( dt, cond );
405        }
406        return ll;
407    }
[488]408
[1064]409    vec samplecond ( const vec &cond );
[488]410
[1064]411    //! }
412    //! \endcode
413    //!@}
[1068]414   
415    /*! Create object from the following structure
416    \code
417    class = 'mmix';
418    pdfs = { list of components bdm::pdf };          % list of pdf offsprings, bdm::pdf::from_setting
419    --- optional fields ---
420    weights = [...];                                 % weights of pdfs in the mixture
421    --- inherited fields ---
422    bdm::pdf::from_setting
423    \endcode
424    \endcode
425    If the optional fields are not given, they will be filled as follows:
426    \code
427    weights = 1/n * [1,1,1,...];
428    \endcode
429    */
[1064]430    void from_setting ( const Setting &set );
[1066]431    void     to_setting  (Setting  &set) const;
[1064]432    virtual void validate();
[488]433};
[737]434SHAREDPTR ( mmix );
[711]435UIREGISTER ( mmix );
[488]436
[979]437
[1066]438//! \brief Base class for all BM running as parallel update of internal BMs
[979]439class ProdBMBase : public BM {
[1064]440protected :
441    Array<vec_from_vec> bm_yt;
442    Array<vec_from_2vec> bm_cond;
[1068]443
444    //! \brief Internal class
[1064]445    class eprod_bm : public eprod_base {
446        ProdBMBase & pb;
447    public :
448        eprod_bm(ProdBMBase &pb0): pb(pb0) {}
449        const epdf* factor(int i ) const {
450            return &(pb.bm(i)->posterior());
451        }
452        const int no_factors() const {
453            return pb.no_bms();
454        }
455    } est;
456public:
457    ProdBMBase():est(*this) {}
458    virtual BM* bm(int i) NOT_IMPLEMENTED(NULL);
459    virtual int no_bms() const {
460        return 0;
461    }
462    const epdf& posterior() const {
463        return est;
464    }
465    void set_prior(const epdf *pri) {
466        const eprod_base* ep=dynamic_cast<const eprod_base*>(pri);
467        if (ep) {
468            bdm_assert(ep->no_factors()!=no_bms() , "Given prior has "+ num2str(ep->no_factors()) + " while this ProdBM has "+
469                       num2str(no_bms()) + "BMs");
470            for (int i=0; i<no_bms(); i++) {
471                bm(i)->set_prior(ep->factor(i));
472            }
473        }
474    }
475
476    void validate() {
477        est.validate();
478        BM::validate();
479        // set links
480        bm_yt.set_length(no_bms());
481        bm_cond.set_length(no_bms());
482
483        //
484
485        for (int i=0; i<no_bms(); i++) {
486            yrv.add(bm(i)->_yrv());
487            rvc.add(bm(i)->_rvc());
488        }
489        rvc=rvc.subt(yrv);
490
491        dimy = yrv._dsize();
492        dimc = rvc._dsize();
493
494        for (int i=0; i<no_bms(); i++) {
495            bm_yt(i).connect(bm(i)->_yrv(), yrv);
496            bm_cond(i).connect(bm(i)->_rvc(), yrv, rvc);
497        }
498    }
499    void bayes(const vec &dt, const vec &cond) {
500        ll=0;
501        for(int i=0; i<no_bms(); i++) {
502            bm_yt(i).update(dt);
503            bm_cond(i).update(dt,cond);
504            bm(i)->bayes(bm_yt(i), bm_cond(i));
505        }
506    }
507
[979]508};
509
[1064]510class ProdBM: public ProdBMBase {
511protected:
512    Array<shared_ptr<BM> > BMs;
513public:
514    virtual BM* bm(int i) {
515        return BMs(i).get();
516    }
517    virtual int no_bms() const {
518        return BMs.length();
519    }
520    void from_setting(const Setting &set) {
521        BM::from_setting(set);
522        UI::get(BMs,set,"BMs");
523    }
524    void to_setting(Setting &set) const {
525        BM::to_setting(set);
526        UI::save(BMs,set,"BMs");
527    }
528
[979]529};
530UIREGISTER(ProdBM);
531
[1064]532//! \brief class for on-line model selection
533class ModelComparator: public ProdBM {
534protected:
535    multiBM weights;
536public:
537    void bayes(const vec &yt, const vec &cond) {
538        vec w_nn(no_bms());
539        for (int i=0; i<no_bms(); i++) {
540            bm(i)->bayes(yt,cond);
541            w_nn(i) += bm(i)->_ll();
542        }
543        vec w=exp(w_nn-max(w_nn));
544        weights.bayes(w/sum(w));
545    }
546    void validate() {
547        ProdBM::validate();
548        weights.validate();
549    }
550    void from_setting(const Setting& set) {
551        ProdBM::from_setting(set);
552        UI::get(weights.frg, set, "frg",UI::optional);
553    }
554    void to_setting(Setting& set) const {
555        ProdBM::to_setting(set);
556        UI::save(weights.frg, set, "frg");
557    }
[1020]558};
559
[254]560}
[107]561#endif //MX_H
Note: See TracBrowser for help on using the browser.