Work at SourceForge, help us to make it a better place! We have an immediate need for a Support Technician in our San Francisco or Denver office.

Close

Commit [669123] default Maximize Restore History

Move RScalarDist and associated D-P-Q-Functions into the main library.

martyn martyn 2014-01-10

1 2 3 > >> (Page 1 of 3)
removed src/modules/bugs/functions/DPQFunction.h
removed src/modules/bugs/functions/PFunction.cc
removed src/modules/bugs/functions/PFunction.h
removed src/modules/bugs/functions/QFunction.cc
removed src/modules/bugs/distributions/DBetaBin.cc
changed src/lib/function/Makefile.am
changed src/lib/distribution/Makefile.am
changed src/lib/Module.cc
changed src/modules/bugs/functions/testbugsfun.cc
changed src/modules/bugs/functions/Makefile.am
changed src/modules/bugs/distributions/DExp.h
changed src/modules/bugs/distributions/DLogis.h
changed src/modules/bugs/distributions/DBeta.h
changed src/modules/bugs/distributions/DGenGamma.h
changed src/modules/bugs/distributions/DWeib.h
copied src/modules/bugs/functions/DFunction.cc -> src/lib/function/DFunction.cc
copied src/modules/bugs/functions/QFunction.h -> src/lib/function/DPQFunction.cc
copied src/modules/bugs/functions/DFunction.h -> src/lib/function/QFunction.cc
copied src/modules/bugs/functions/DPQFunction.cc -> src/lib/function/PFunction.cc
copied src/modules/bugs/distributions/RScalarDist.h -> src/lib/distribution/RScalarDist.cc
src/lib/function/Makefile.am Diff Switch to side-by-side view
Loading...
src/lib/distribution/Makefile.am Diff Switch to side-by-side view
Loading...
src/lib/Module.cc Diff Switch to side-by-side view
Loading...
src/modules/bugs/functions/testbugsfun.cc Diff Switch to side-by-side view
Loading...
src/modules/bugs/functions/Makefile.am Diff Switch to side-by-side view
Loading...
src/modules/bugs/distributions/DExp.h Diff Switch to side-by-side view
Loading...
src/modules/bugs/distributions/DLogis.h Diff Switch to side-by-side view
Loading...
src/modules/bugs/distributions/DBeta.h Diff Switch to side-by-side view
Loading...
src/modules/bugs/distributions/DGenGamma.h Diff Switch to side-by-side view
Loading...
src/modules/bugs/distributions/DWeib.h Diff Switch to side-by-side view
Loading...
src/modules/bugs/functions/DFunction.cc to src/lib/function/DFunction.cc
--- a/src/modules/bugs/functions/DFunction.cc
+++ b/src/lib/function/DFunction.cc
@@ -1,12 +1,11 @@
 #include <config.h>
-#include "RScalarDist.h"
-#include "DFunction.h"
+#include <distribution/RScalarDist.h>
+#include <function/DFunction.h>
 
 using std::vector;
 using std::string;
 
 namespace jags {
-namespace bugs {
 
     DFunction::DFunction(RScalarDist const *dist)
 	: DPQFunction(dist->name(), dist)
@@ -35,4 +34,4 @@
 	return checkArgs(args);
     }
 
-}}
+}
src/modules/bugs/functions/QFunction.h to src/lib/function/DPQFunction.cc
--- a/src/modules/bugs/functions/QFunction.h
+++ b/src/lib/function/DPQFunction.cc
@@ -1,19 +1,29 @@
-#ifndef Q_FUNCTION_H_
-#define Q_FUNCTION_H_
+#include <config.h>
+#include <function/DPQFunction.h>
+#include <distribution/RScalarDist.h>
 
-#include "DPQFunction.h"
+using std::vector;
+using std::string;
 
 namespace jags {
-namespace bugs {
 
-    class QFunction : public DPQFunction
+    DPQFunction::DPQFunction(string const &name, RScalarDist const *dist) 
+	: ScalarFunction(name, dist->npar() + 1), _dist(dist)
+    {}
+    
+    RScalarDist const *DPQFunction::dist() const
     {
-    public:
-	QFunction(RScalarDist const *dist);
-	bool checkParameterValue(std::vector<double const *> const &args) const;
-	double evaluate(std::vector <double const *> const &args) const;
-    };
+	return _dist;
+    }
+    
+    bool DPQFunction::checkArgs(vector<double const *> const &args) const
+    {
+	vector<double const *> param(_dist->npar());
+	for (unsigned int i = 0; i < param.size(); ++i) {
+	    param[i] = args[i+1];
+	}
+	
+	return _dist->checkParameterValue(param);
+    }
 
-}}
-
-#endif /* Q_FUNCTION_H_ */
+}
src/modules/bugs/functions/DFunction.h to src/lib/function/QFunction.cc
--- a/src/modules/bugs/functions/DFunction.h
+++ b/src/lib/function/QFunction.cc
@@ -1,19 +1,30 @@
-#ifndef D_FUNCTION_H_
-#define D_FUNCTION_H_
+#include <config.h>
+#include <function/QFunction.h>
+#include <distribution/RScalarDist.h>
 
-#include "DPQFunction.h"
+using std::vector;
+using std::string;
 
 namespace jags {
-namespace bugs {
 
-    class DFunction : public DPQFunction
+    QFunction::QFunction(RScalarDist const *dist)
+	: DPQFunction(string("q") + dist->name().substr(1), dist)
+    {}
+    
+    double QFunction::evaluate(vector<double const *> const &args) const
     {
-    public:
-	DFunction(RScalarDist const *dist);
-	bool checkParameterValue(std::vector<double const *> const &args) const;
-	double evaluate(std::vector <double const *> const &args) const;
-    };
+	double x = *args[0];
+	vector<double const *> param(args.size() - 1);
+	for (unsigned int i = 1; i < args.size(); ++i) {
+	    param[i-1] = args[i];
+	}
+	
+	return dist()->q(x, param, true, false);
+    }
 
-}}
+    bool QFunction::checkParameterValue(vector<double const*> const &args) const
+    {
+	return checkArgs(args);
+    }
 
-#endif /* D_FUNCTION_H_ */
+}
src/modules/bugs/functions/DPQFunction.cc to src/lib/function/PFunction.cc
--- a/src/modules/bugs/functions/DPQFunction.cc
+++ b/src/lib/function/PFunction.cc
@@ -1,30 +1,37 @@
 #include <config.h>
-#include "DPQFunction.h"
-#include "RScalarDist.h"
+#include <function/PFunction.h>
+#include <distribution/RScalarDist.h>
 
 using std::vector;
 using std::string;
 
 namespace jags {
-namespace bugs {
 
-    DPQFunction::DPQFunction(string const &name, RScalarDist const *dist) 
-	: ScalarFunction(name, dist->npar() + 1), _dist(dist)
+    PFunction::PFunction(RScalarDist const *dist)
+	: DPQFunction(string("p") + dist->name().substr(1), dist)
     {}
     
-    RScalarDist const *DPQFunction::dist() const
+    double PFunction::evaluate(vector<double const *> const &args) const
     {
-	return _dist;
-    }
-    
-    bool DPQFunction::checkArgs(vector<double const *> const &args) const
-    {
-	vector<double const *> param(_dist->npar());
-	for (unsigned int i = 0; i < param.size(); ++i) {
-	    param[i] = args[i+1];
+	double x = *args[0];
+	vector<double const *> param(args.size() - 1);
+	for (unsigned int i = 1; i < args.size(); ++i) {
+	    param[i-1] = args[i];
 	}
 	
-	return _dist->checkParameterValue(param);
+	return dist()->p(x, param, true, false);
     }
 
-}}
+    bool 
+    PFunction::checkParameterValue(vector<double const *> const &args) const
+    {
+	if (dist()->discrete()) {
+	    double x = *args[0];
+	    if (x != static_cast<int>(x))
+		return false;
+	}
+
+	return checkArgs(args);
+    }
+
+}
src/modules/bugs/distributions/RScalarDist.h to src/lib/distribution/RScalarDist.cc
--- a/src/modules/bugs/distributions/RScalarDist.h
+++ b/src/lib/distribution/RScalarDist.cc
@@ -1,125 +1,186 @@
-#ifndef R_SCALAR_DIST_H_
-#define R_SCALAR_DIST_H_
+#include <config.h>
+#include <distribution/RScalarDist.h>
+#include <rng/RNG.h>
+#include <util/nainf.h>
+#include <util/dim.h>
 
-#include <distribution/ScalarDist.h>
+#include <cmath>
+#include <algorithm>
+
+using std::string;
+using std::vector;
+using std::log;
+using std::min;
+using std::max;
 
 namespace jags {
 
-struct RNG;
+double RScalarDist::calPlower(double lower, 
+			      vector<double const*> const &parameters) const
+{
+    //P(X < lower)
+    if (_discrete) lower -= 1;
+    return p(lower, parameters, true, false);
+}
 
-namespace bugs {
+double RScalarDist::calPupper(double upper,
+			      vector<double const*> const &parameters) const
+{
+    //P(X <= upper)
+    return p(upper, parameters, true, false);
+}
 
-/**
- * @short Scalar Distribution using R math library infrastructure.
- *
- * A subclass of RScalarDist has to implement the d,p,q, and r virtual
- * member functions. These are based on the d-p-q-r functions provided
- * by libRmath.
- *
- * The JAGS versions of most (but not all) scalar distributions extend
- * the distribution families in libRmath by allowing the distribution
- * to be bounded.
- */
-class RScalarDist : public ScalarDist
+
+RScalarDist::RScalarDist(string const &name, unsigned int npar, 
+			 Support support, bool discrete)
+  
+    : ScalarDist(name, npar, support),  _support(support), _discrete(discrete),
+      _npar(npar)
 {
-    const Support _support;
-    const bool _discrete;
-    unsigned int _npar;
-    double calPlower(double, std::vector<double const *> const &) const;
-    double calPupper(double, std::vector<double const *> const &) const;
-public:
-    /**
-     * Constructor
-     *
-     * @param name BUGS language name of distribution
-     *
-     * @param npar Number of parameters, excluding upper and lower bound
-     *
-     * @param support Support of distribution
-     *
-     * @param discrete Boolean flag indicating whether the distribution is
-     *        discrete valued.
-     */
-    RScalarDist(std::string const &name, unsigned int npar, Support support,
-		bool discrete=false);
-    double logDensity(double x, PDFType type,
-		      std::vector<double const *> const &parameters,
-		      double const *lower, double const *upper) const;
-    double randomSample(std::vector<double const *> const &parameters,
-			double const *lower, double const *upper,
-			RNG *rng) const;
-    /**
-     * Returns the median. Note that this function can be overloaded
-     * by a subclass if necessary.
-     */
-    double typicalValue(std::vector<double const *> const &parameters,
-			double const *lower, double const *upper) const;
-    /**
-     * Density function, ignoring bounds
-     * @param x value at which to evaluate the density
-     * @param type Type of density calculation required.
-     * @param parameters Array of parameters
-     * @param give_log Indicates whether to return log density. 
-     */
-    virtual double d(double x, PDFType type,
-		     std::vector<double const *> const &parameters, 
-		     bool give_log) const = 0;
-    /**
-     * Distribution function, ignoring bounds
-     * @param x quantile at which to evaluate the distribution function
-     * @param parameters Array of parameters
-     * @param lower If true, return value is P[X <= x]. Otherwise
-     * P[X > x]
-     * @param give_log Indicates whether to return log probabability
-     */
-    virtual double p(double x, std::vector<double const *> const &parameters, 
-		     bool lower, bool give_log) const = 0;
-    /**
-     * Quantile function, ignoring bounds
-     * @param p probability for which to evaluate quantile
-     * @param parameters Array of parameters
-     * @param log_p Indicates whether p is given as log(p). 
-     */
-    virtual double q(double p, std::vector<double const *> const &parameters, 
-		     bool lower, bool log_p) const = 0;
-    /**
-     * Random number generation, ignoring bounds
-     * @param parameters Array of parameters
-     */
-    virtual double 
-	r(std::vector<double const *> const &parameters, RNG *rng) const = 0;
-    /**
-     * All RScalarDist distributions can be bounded
-     */
-    bool canBound() const;
-    /**
-     * RScalarDist distributions are defined to have support on the integers
-     * or on the real line by the constructor
-     */
-    bool isDiscreteValued(std::vector<bool> const &mask) const;
-    /**
-     * Alternative function for determining whether the distribution is
-     * discrete-valued.
-     */
-    bool discrete() const;
-    /**
-     * Returns the number of parameters of the distribution
-     */
-    unsigned int npar() const;
-};
+}
 
-    /**
-     * Convenience function that calculates x * log(0) as the limit of
-     * x * log(p) as p tends to zero.  This is required for calculation
-     * of some density functions.
-     *
-     * @param x coefficient of log(0)
-     * @param give_log logical flag. If true then the limit of x*log(p)
-     * is returned, otherwise the limit of p^x.
-     */
-    double xlog0(double x, bool give_log);
+double 
+RScalarDist::typicalValue(vector<double const *> const &parameters,
+			  double const *lower, double const *upper) const
+{
+    double llimit = l(parameters), ulimit = u(parameters);
+    double plower = 0, pupper = 1;
+    
+    if (lower) {
+	llimit = max(llimit, *lower);
+	plower = calPlower(llimit, parameters);
+    }
 
-}}
+    if (upper) {
+	ulimit = min(ulimit, *upper);
+	pupper = calPupper(ulimit, parameters);
+    }
+    
+    double pmed = (plower + pupper)/2;
+    double med = q(pmed, parameters, true, false);	
 
-#endif /* SCALAR_DIST_RMATH_H_ */
+    //Calculate the log densities
+    double dllimit = d(llimit, PDF_FULL, parameters, true);
+    double dulimit = d(ulimit, PDF_FULL, parameters, true);
+    double dmed = d(med, PDF_FULL, parameters, true);
 
+    //Pick the median if it has the highest density, otherwise pick
+    //a point near to (but not on) the boundary
+    if (dmed >= dllimit && dmed >= dulimit) {
+	return med;
+    }
+    else if (dulimit > dllimit) {
+	return q(0.1 * plower + 0.9 * pupper, parameters, true, false);
+    }
+    else {
+	return q(0.9 * plower + 0.1 * pupper, parameters, true, false);
+    }
+}
+
+double 
+RScalarDist::logDensity(double x, PDFType type,
+			vector<double const *> const &parameters,
+			double const *lower, double const *upper) const
+{
+    if (lower && x < *lower)
+	return JAGS_NEGINF;
+    if (upper && x > *upper)
+	return JAGS_NEGINF;
+    if (upper && lower && *upper < *lower)
+	return JAGS_NEGINF;
+    
+    double loglik =  d(x, type, parameters, true);
+
+    if (type != PDF_PRIOR && (lower || upper)) {
+	//Normalize truncated distributions
+
+	double ll = l(parameters);
+	if (lower && *lower < ll) ll = *lower;
+	if (_discrete) ll -= 1; //Adjustment for discrete valued distributions
+
+	/* In theory, we just have to subtract log[P(lower <= X <=
+           upper)] from the log likelihood. But we need to work around
+           numerical problems. */
+
+	bool have_lower = lower && p(ll, parameters, true, false) > 0;
+	bool have_upper = upper && p(*upper, parameters, false, false) > 0;
+
+	if (have_lower && have_upper) {
+	    if (p(ll, parameters, false, false) < 0.5) {
+		//Use upper tail
+		loglik -= log(p(ll, parameters, false, false) -
+			      p(*upper, parameters, false, false));
+	    }
+	    else {
+		//Use lower tail
+		loglik -= log(p(*upper, parameters, true, false) - 
+			      p(ll, parameters, true, false));
+	    }
+	}
+	else if (have_lower) {
+	    loglik -= p(ll, parameters, false, true);
+	}
+	else if (have_upper) {
+	    loglik -= p(*upper, parameters, true, true);
+	}
+    }
+
+    return loglik;
+}
+
+
+double 
+RScalarDist::randomSample(vector<double const *> const &parameters,
+			  double const *lower, double const *upper,
+			  RNG *rng) const
+{
+    if (!lower && !upper) {
+	return r(parameters, rng);
+    }
+
+    double plower = lower ? calPlower(*lower, parameters) : 0;
+    double pupper = upper ? calPupper(*upper, parameters) : 1;
+
+    if (pupper - plower > 0.25) {
+	//Rejection sampling if expected number of samples is 4 or less
+	while (true) {
+	    double y = r(parameters, rng);
+	    if (lower && y < *lower) continue;
+	    if (upper && y > *upper) continue;
+	    return y;
+	}
+    }
+
+    //Inversion
+    //FIXME: We probably need to take care of tail behaviour here
+    double u = plower + rng->uniform() * (pupper - plower);
+    return q(u, parameters, true, false);
+}
+
+bool RScalarDist::canBound() const
+{
+    return true;
+}
+
+bool RScalarDist::isDiscreteValued(vector<bool> const &mask) const
+{
+    return _discrete;
+}
+
+bool RScalarDist::discrete() const
+{
+    return _discrete;
+}
+
+unsigned int RScalarDist::npar() const
+{
+    return _npar;
+}
+
+    double xlog0(double x, bool give_log) {
+	if (x < 0) return JAGS_POSINF;
+	else if (x > 0) return give_log ? JAGS_NEGINF : 0;
+	else return give_log ? 0 : 1;
+    }
+
+}
1 2 3 > >> (Page 1 of 3)