src/prior.h

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

#include <cmath>
#include "rng.h"
#include "special.h"

class PsiPrior
{
        private:
                PsiRandom rng;
        public:
                virtual double pdf ( double x ) const { return 1.;}    
                virtual double dpdf ( double x ) { return 0.; }  
                virtual double rand ( void ) { return rng.draw(); } 
                virtual PsiPrior * clone ( void ) const { throw NotImplementedError(); }
                virtual double mean ( void ) const { return 0; } 
                virtual double std  ( void ) const { return 1e5; } 
                virtual void shrink ( double xmin, double xmax ) { throw NotImplementedError(); } 
                virtual int get_code(void) const { throw NotImplementedError(); } 
                virtual double cdf ( double x ) const { throw NotImplementedError(); } 
                virtual double getprm ( unsigned int prm ) const { throw NotImplementedError(); }
                virtual double ppf ( double p, double start=NULL ) const { throw NotImplementedError(); }
};

class UniformPrior : public PsiPrior
{
        private:
                double lower;
                double upper;
                double height;
                UniformRandom rng;
        public:
                UniformPrior ( double low, double high ) : lower(low), upper(high), rng ( low, high ) { height = 1./(high-low); } 
                UniformPrior ( const UniformPrior& original ) : lower(original.lower),
                                                        upper(original.upper),
                                                        height(original.height),
                                                        rng(original.rng) {} 
                double pdf ( double x ) const { return ( x>lower && x<upper ? height : 0 ); }                      
                double dpdf ( double x ) { return ( x!=lower && x!=upper ? 0 : (x==lower ? 1e20 : -1e20 ));} 
                double rand ( void ) { return rng.draw(); }                                                 
        PsiPrior * clone ( void ) const { return new UniformPrior(*this); }
                double mean ( void ) const { return 0.5*(lower+upper); }  
                double std  ( void ) const { return sqrt((upper-lower)*(upper-lower)/12); }
                void shrink ( double xmin, double xmax ) {}         
                int get_code(void) const { return 0; } 
                double cdf ( double x ) const { return ( x<lower ? 0 : (x>upper ? 1 : (x-lower)/(upper-lower) ) ); }
                double getprm ( unsigned int prm ) const { return (prm==0 ? lower : upper ); }
                double ppf ( double p, double start=NULL ) const { return ( p>1 ? upper : (p<0 ? lower : p*(upper-lower)+lower)); }
};

class GaussPrior : public PsiPrior
{
        private:
                double mu;
                double sg;
                double normalization;
                double var;
                double twovar;
                GaussRandom rng;
        public:
                GaussPrior ( double mean, double sd ) : mu(mean), sg(sd), var(sg*sg), twovar(2*sg*sg), rng(mean,sd) { normalization = 1./(sqrt(2*M_PI)*sg); }        
                GaussPrior ( const GaussPrior& original ) : mu(original.mu),
                                                    sg(original.sg),
                                                    normalization(original.normalization),
                                                    var(original.var),
                                                    twovar(original.twovar),
                                                    rng(original.rng) {} 
                double pdf ( double x ) const { return normalization * exp ( - (x-mu)*(x-mu)/twovar ); }                                              
                double dpdf ( double x ) { return - x * pdf ( x ) / var; }                                                                      
                double rand ( void ) {return rng.draw(); }
        PsiPrior * clone ( void ) const { return new GaussPrior(*this); }
                double mean ( void ) const { return mu; } 
                double std  ( void ) const { return sg; } 
                void shrink ( double xmin, double xmax );
                int get_code(void) const { return 1; } 
                double cdf ( double x ) const { return Phi ( (x-mu)/sg ); }
                double getprm ( unsigned int prm ) const { return ( prm==0 ? mu : sg ); }
                double ppf ( double p, double start=NULL ) const;
};

class BetaPrior : public PsiPrior
{
        private:
                double alpha;
                double beta;
                double normalization;
                BetaRandom rng;
                double mode;
        public:
                BetaPrior ( double al, double bt ) : alpha(al), beta(bt), normalization(betaf(al,bt)), rng (alpha, beta) {
            mode = (al-1)/(al+bt-2);
            mode = pdf(mode); }                      
        BetaPrior ( const BetaPrior& original) : alpha(original.alpha),
                                                 beta(original.beta),
                                                 normalization(original.normalization),
                                                 rng(original.rng),
                                                 mode(original.mode) {} 
                double pdf ( double x ) const { return (x<0||x>1 ? 0 : pow(x,alpha-1)*pow(1-x,beta-1)/normalization); }             
                double dpdf ( double x ) { return (x<0||x>1 ? 0 : ((alpha-1)*pow(x,alpha-2)*pow(1-x,beta-1) + (beta-1)*pow(1-x,beta-2)*pow(x,alpha-1))/normalization); }      
                double rand ( void ) {return rng.draw();};                                                                                         
        PsiPrior * clone ( void ) const { return new BetaPrior(*this); }
                double mean ( void ) const { return alpha/(alpha+beta); }
                double std  ( void ) const { return sqrt ( alpha*beta/((alpha+beta)*(alpha+beta)*(alpha+beta+1)) ); }
                void shrink ( double xmin, double xmax );
                int get_code(void) const { return 2; } 
                double cdf ( double x ) const { return (x<0 ? 0 : (x>1 ? 1 : betainc ( x, alpha, beta ))); }
                double getprm ( unsigned int prm ) const { return ( prm==0 ? alpha : beta ); }
                double ppf ( double p, double start=NULL ) const;
};

class GammaPrior : public PsiPrior
{
        private:
                double k;
                double theta;
                double normalization;
                GammaRandom rng;
        public:
                GammaPrior ( double shape, double scale ) : k(shape), theta(scale), rng(shape, scale) { normalization = pow(scale,shape)*exp(gammaln(shape));}                         
        GammaPrior ( const GammaPrior& original ) : k(original.k),
                                                    theta(original.theta),
                                                    normalization(original.normalization),
                                                    rng(original.rng) {} 
                virtual double pdf ( double x ) const { return (x>0 ? pow(x,k-1)*exp(-x/theta)/normalization : 0 );}                                                             
                virtual double dpdf ( double x ) { return (x>0 ? ( (k-1)*pow(x,k-2)*exp(-x/theta)-pow(x,k-1)*exp(-x/theta)/theta)/normalization : 0 ); }                   
                virtual double rand ( void ) {return rng.draw(); };
        PsiPrior * clone ( void ) const { return new GammaPrior(*this); }
                virtual double mean ( void ) const { return k*theta; }
                double std  ( void ) const { return sqrt ( k*theta*theta ); }
                void shrink ( double xmin, double xmax );
                virtual int get_code(void) const { return 3; } 
                virtual double cdf ( double x ) const { return ( x<0 ? 0 : gammainc ( k, x/theta ) / exp ( gammaln ( k ) ) ); }
                double getprm ( unsigned int prm ) const { return ( prm==0 ? k : theta ); }
                virtual double ppf ( double p, double start=NULL ) const;
};

class nGammaPrior : public GammaPrior
{
        public:
                nGammaPrior ( double shape, double scale ) : GammaPrior(shape,scale) {}
        nGammaPrior ( const nGammaPrior& original ) : GammaPrior(original) {} 
                double pdf ( double x ) const { return GammaPrior::pdf ( -x ); }
                double dpdf ( double x ) { return -GammaPrior::dpdf ( -x ); }
                double rand ( void ) { return -GammaPrior::rand(); }
        PsiPrior * clone ( void ) const { return new nGammaPrior(*this); }
                double mean ( void ) const { return -GammaPrior::mean(); }
                void shrink ( double xmin, double xmax );
                int get_code(void) const { return 4; } 
                double cdf ( double x ) const { return ( x>0 ? 1 : 1-GammaPrior::cdf ( -x ) ); }
                double ppf ( double p, double start=NULL ) const { return - GammaPrior::ppf ( 1-p ); }
};

class invGammaPrior : public PsiPrior
{
        private:
                double alpha;
                double beta;
                double normalization;
                GammaRandom rng;
        public:
                invGammaPrior ( double shape, double scale ) : alpha ( shape ), beta ( scale ), rng ( shape, 1./scale ) { normalization = pow(alpha,beta)/exp(gammaln(shape));} 
                invGammaPrior ( const invGammaPrior& original ) : alpha(original.alpha),
                                        beta(original.beta),
                                        normalization ( original.normalization ),
                                        rng(original.rng) {} 
                virtual double pdf ( double x ) { return ( x>0 ? pow ( x, -alpha-1 ) * exp ( -beta/x ) * normalization : 0 ); }
                virtual double dpdf ( double x ) { return (x>0 ? ( (-alpha-1)*pow(x,-alpha-2) * exp ( -beta/x ) + pow(x,-alpha-1) * exp ( -beta/x ) * beta / (x*x) ) * normalization : 0 ); }
                virtual double rand ( void ) { return 1./rng.draw(); }
                PsiPrior * clone ( void ) const { return new invGammaPrior(*this); }
                virtual double mean ( void ) const { return beta/(alpha-1); }
                double std ( void ) const { return ( alpha>2 ? beta / ( (alpha-1)*sqrt(alpha-2) ) : 1e5 ); }
                virtual void shrink ( double xmin, double xmax ) {} 
                virtual int get_code ( void ) const { return 5; } 
};

class ninvGammaPrior : public invGammaPrior
{
        public:
                ninvGammaPrior ( double shape, double scale ) : invGammaPrior ( shape, scale ) {}
                ninvGammaPrior ( const ninvGammaPrior& original ) : invGammaPrior ( original ) {}
                double pdf ( double x ) { return invGammaPrior::pdf ( -x ); }
                double dpdf ( double x ) { return -invGammaPrior::dpdf ( -x ); }
                double rand ( void ) { return -invGammaPrior::rand(); }
                PsiPrior * clone ( void ) const { return new ninvGammaPrior ( *this ); }
                double mean ( void ) const { return -invGammaPrior::mean(); }
                void shrink ( double xmin, double xmax ) { invGammaPrior::shrink ( -xmax, -xmin ); }
                int get_code ( void ) const { return 6; } 
};
#endif