src/mcmc.h

/*
 *   See COPYING file distributed along with the psignifit package for
 *   the copyright and license terms
 */
#ifndef MCMC_H
#define MCMC_H

#include <vector>
#include "psychometric.h"
#include "rng.h"
#include "mclist.h"
#include "getstart.h"

class PsiSampler
{
        private:
                const PsiPsychometric * model;
                const PsiData * data;
        public:
                PsiSampler ( const PsiPsychometric * Model, const PsiData * Data ) : model(Model), data(Data) {}
                virtual std::vector<double> draw ( void ) { throw NotImplementedError(); }                     
                virtual void setTheta ( const std::vector<double> theta ) { throw NotImplementedError(); }     
                virtual std::vector<double> getTheta ( void ) { throw NotImplementedError(); }                 
                virtual void setStepSize ( double size, unsigned int param ) { throw NotImplementedError(); }  
                virtual void setStepSize ( const std::vector<double>& sizes ) { throw NotImplementedError(); } 
                virtual double getDeviance ( void ) { throw NotImplementedError(); }                           
                virtual MCMCList sample ( unsigned int N ) { throw NotImplementedError(); }                   
                const PsiPsychometric * getModel() const { return model; }                                     
                const PsiData         * getData()  const { return data;  }                                     
};

class MetropolisHastings : public PsiSampler
{
        private:
                PsiRandom* propose;
                std::vector<double> currenttheta;
                std::vector<double> newtheta;
                std::vector<double> stepwidths;
                double currentdeviance;
                int accept;
        protected:
                double qold;
        public:
                MetropolisHastings (
                        const PsiPsychometric * Model,                                                  
                        const PsiData * Data,                                                           
                        PsiRandom* proposal                                                             
                        );                                                          
                ~MetropolisHastings ( void ) { delete propose; }
                std::vector<double> draw ( void );                                                
                virtual double acceptance_probability ( const std::vector<double>& current_theta, const std::vector<double>& new_theta );
                void setTheta ( const std::vector<double>& prm );                                 
                std::vector<double> getTheta ( void ) { return currenttheta; }                    
                double getDeviance ( void ) { return currentdeviance; }                           
                void setStepSize ( double size, unsigned int param );                             
                void setStepSize ( const std::vector<double>& sizes );                            
                std::vector<double> getStepsize ( void ) { return stepwidths; }                                   
                MCMCList sample ( unsigned int N );                                               
                unsigned int getNparams ( void ) { return newtheta.size(); }                      
                virtual void proposePoint( std::vector<double> &current_theta,
                                                                        std::vector<double> &step_widths,
                                                                        PsiRandom * proposal,
                                                                        std::vector<double> &new_theta);                                          
};

class GenericMetropolis : public MetropolisHastings
{
        private:
                int currentindex;                                                                                                                                 
        public:
                GenericMetropolis (
                        const PsiPsychometric * Model,                                                            
                        const PsiData * Data,                                                                     
                        PsiRandom* proposal                                                                       
                        ): MetropolisHastings ( Model, Data, proposal ),
                           currentindex(0) {}
                void proposePoint( std::vector<double> &current_theta,
                                                        std::vector<double> &step_widths,
                                                        PsiRandom * proposal,
                                                        std::vector<double> &new_theta);                                                          

                void findOptimalStepwidth ( PsiMClist const &pilot );
};

class DefaultMCMC : public MetropolisHastings
{
        private:
                std::vector<PsiPrior*> proposaldistributions;
        public:
                DefaultMCMC ( const PsiPsychometric * Model,                                      
                                const PsiData * Data,                                                     
                                PsiRandom* proposal                                                       
                                );
                ~DefaultMCMC ( void );
                double acceptance_probability (
                                const std::vector<double> &current_theta,
                                const std::vector<double> &new_theta );
                void proposePoint ( std::vector<double> &current_theta,
                                std::vector<double> &stepwidths,
                                PsiRandom * proposal,
                                std::vector<double> &new_theta);                                          
        void set_proposal(unsigned int i, PsiPrior* proposal){
            delete proposaldistributions.at(i);
            proposaldistributions.at(i) = proposal->clone();
        }
};

class HybridMCMC : public PsiSampler
{
        private:
                PsiRandom* proposal;
                std::vector<double> currenttheta;
                std::vector<double> newtheta;
                std::vector<double> momentum;
                double currentH;
                double newH;
                double energy;
                double newenergy;
                std::vector<double> gradient;
                std::vector<double> currentgradient;
                std::vector<double> stepsizes;
                int Nleapfrog;
                int Naccepted;
                void leapfrog ( void );
        public:
                HybridMCMC (
                        const PsiPsychometric * Model,                                                  
                        const PsiData * Data,                                                          
                        int Nleap                                                                     
                        );                                                             
                ~HybridMCMC ( void ) { delete proposal; }
                std::vector<double> draw ( void );                                                
                void setTheta ( const std::vector<double>& prm );                                 
                std::vector<double> getTheta ( void ) { return currenttheta; }                    
                void setStepSize ( double size, unsigned int param );                             
                void setStepSize ( const std::vector<double>& sizes );                            
                double getDeviance ( void );                                                      
                MCMCList sample ( unsigned int N );                                              
};

double ModelEvidence ( const PsiPsychometric* pmf, const PsiData* data );

std::vector<double> OutlierDetection ( const PsiPsychometric* pmf, OutlierModel* outl, const PsiData* data );

#endif