[Rcpp-commits] r2011 - in pkg/Rcpp: . inst/include/Rcpp/stats
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Aug 15 13:08:51 CEST 2010
Author: romain
Date: 2010-08-15 13:08:51 +0200 (Sun, 15 Aug 2010)
New Revision: 2011
Modified:
pkg/Rcpp/TODO
pkg/Rcpp/inst/include/Rcpp/stats/nbinom_mu.h
pkg/Rcpp/inst/include/Rcpp/stats/norm.h
Log:
dpq norm now uses RCPP_DPQ macros
Modified: pkg/Rcpp/TODO
===================================================================
--- pkg/Rcpp/TODO 2010-08-15 10:20:14 UTC (rev 2010)
+++ pkg/Rcpp/TODO 2010-08-15 11:08:51 UTC (rev 2011)
@@ -52,9 +52,7 @@
o not sure rep should be lazy, i.e. rep( x, 4 ) fetches x[i] 4 times,
maybe we should use LazyVector like in outer to somehow cache the
result when it is judged expensive to calculate
-
- o convert dpq norm and dpq unif to use RCPP_DPQ
-
+
o crossprod
o SUGAR_MATH: is there overhead in having the function pointer, should
Modified: pkg/Rcpp/inst/include/Rcpp/stats/nbinom_mu.h
===================================================================
--- pkg/Rcpp/inst/include/Rcpp/stats/nbinom_mu.h 2010-08-15 10:20:14 UTC (rev 2010)
+++ pkg/Rcpp/inst/include/Rcpp/stats/nbinom_mu.h 2010-08-15 11:08:51 UTC (rev 2011)
@@ -22,7 +22,7 @@
#ifndef Rcpp__stats__nbinom_mu_h
#define Rcpp__stats__nbinom_mu_h
-RCPP_DPQ_2(nbinom_mu,::Rf_dnbinom_mu, ::Rf_pnbinom_mu, ::Rf_qnbinom_mu )
+RCPP_DPQ_2(nbinom_mu,::dnbinom_mu, ::pnbinom_mu, ::qnbinom_mu )
#endif
Modified: pkg/Rcpp/inst/include/Rcpp/stats/norm.h
===================================================================
--- pkg/Rcpp/inst/include/Rcpp/stats/norm.h 2010-08-15 10:20:14 UTC (rev 2010)
+++ pkg/Rcpp/inst/include/Rcpp/stats/norm.h 2010-08-15 11:08:51 UTC (rev 2011)
@@ -23,93 +23,88 @@
#define Rcpp__stats__norm_h
namespace Rcpp {
-
namespace stats {
+inline double dnorm_1(double x, double mu /*, double sigma [=1.0]*/ , int give_log) {
+#ifdef IEEE_754
+ if (ISNAN(x) || ISNAN(mu) )
+ return x + mu + 1.0 ;
+#endif
+ if(!R_FINITE(x) && mu == x) return ML_NAN;/* x-mu is NaN */
+ x = (x - mu) ;
+ if(!R_FINITE(x)) return R_D__0;
+ return (give_log ?
+ -(M_LN_SQRT_2PI + 0.5 * x * x ) :
+ M_1_SQRT_2PI * ::exp(-0.5 * x * x) );
+ /* M_1_SQRT_2PI = 1 / sqrt(2 * pi) */
+}
- template <bool NA, typename T>
- class DNorm : public Rcpp::VectorBase< REALSXP, NA, DNorm<NA,T> >{
- public:
- typedef typename Rcpp::VectorBase<REALSXP,NA,T> VEC_TYPE;
-
- DNorm( const VEC_TYPE& vec_, double mu_, double sigma_, bool log_ = false ) :
- vec(vec_), mu(mu_), sigma(sigma_), log(log_) {}
-
- inline double operator[]( int i) const {
- return ::Rf_dnorm4( vec[i], mu, sigma, log );
- }
-
- inline int size() const { return vec.size(); }
-
- private:
- const VEC_TYPE& vec;
- double mu, sigma;
- int log;
-
- };
+inline double dnorm_0(double x /*, double mu [=0.0], double sigma [=1.0]*/ , int give_log) {
+#ifdef IEEE_754
+ if (ISNAN(x) )
+ return x + 1.0 ;
+#endif
+ if(!R_FINITE(x)) return R_D__0;
+ return (give_log ?
+ -(M_LN_SQRT_2PI + 0.5 * x * x ) :
+ M_1_SQRT_2PI * ::exp(-0.5 * x * x) );
+ /* M_1_SQRT_2PI = 1 / sqrt(2 * pi) */
+}
- template <bool NA, typename T>
- class PNorm : public Rcpp::VectorBase< REALSXP, NA, PNorm<NA,T> >{
- public:
- typedef typename Rcpp::VectorBase<REALSXP,NA,T> VEC_TYPE;
-
- PNorm( const VEC_TYPE& vec_, double mu_, double sigma_,
- bool lower_tail = true, bool log_ = false ) :
- vec(vec_), mu(mu_), sigma(sigma_), lower(lower_tail), log(log_) {}
-
- inline double operator[]( int i) const {
- return ::Rf_pnorm5( vec[i], mu, sigma, lower, log );
- }
-
- inline int size() const { return vec.size(); }
-
- private:
- const VEC_TYPE& vec;
- double mu, sigma;
- int lower, log;
-
- };
+inline double pnorm_1(double x, double mu /*, double sigma [=1.]*/ , int lower_tail, int log_p){
+ double p, cp;
- template <bool NA, typename T>
- class QNorm : public Rcpp::VectorBase< REALSXP, NA, QNorm<NA,T> >{
- public:
- typedef typename Rcpp::VectorBase<REALSXP,NA,T> VEC_TYPE;
-
- QNorm( const VEC_TYPE& vec_, double mu_, double sigma_,
- bool lower_tail = true, bool log_ = false ) :
- vec(vec_), mu(mu_), sigma(sigma_), lower(lower_tail), log(log_) {}
-
- inline double operator[]( int i) const {
- return ::Rf_qnorm5( vec[i], mu, sigma, lower, log );
- }
-
- inline int size() const { return vec.size(); }
-
- private:
- const VEC_TYPE& vec;
- double mu, sigma;
- int lower, log;
-
- };
-
-} // stats
+ /* Note: The structure of these checks has been carefully thought through.
+ * For example, if x == mu and sigma == 0, we get the correct answer 1.
+ */
+#ifdef IEEE_754
+ if(ISNAN(x) || ISNAN(mu) )
+ return x + mu + 1.0 ;
+#endif
+ if(!R_FINITE(x) && mu == x) return ML_NAN;/* x-mu is NaN */
+ p = (x - mu) ;
+ if(!R_FINITE(p))
+ return (x < mu) ? R_DT_0 : R_DT_1;
+ x = p;
-template <bool NA, typename T>
-inline stats::DNorm<NA,T> dnorm( const Rcpp::VectorBase<REALSXP,NA,T>& x, double mu, double sigma, bool log = false ) {
- return stats::DNorm<NA,T>( x, mu, sigma, log );
+ ::Rf_pnorm_both(x, &p, &cp, (lower_tail ? 0 : 1), log_p);
+
+ return(lower_tail ? p : cp);
}
-template <bool NA, typename T>
-inline stats::PNorm<NA,T> pnorm( const Rcpp::VectorBase<REALSXP,NA,T>& x, double mu, double sigma, bool lower = true, bool log = false ) {
- return stats::PNorm<NA,T>( x, mu, sigma, lower, log );
+inline double pnorm_0(double x /*, double mu [=0.] , double sigma [=1.]*/ , int lower_tail, int log_p) {
+ double p, cp;
+
+ /* Note: The structure of these checks has been carefully thought through.
+ * For example, if x == mu and sigma == 0, we get the correct answer 1.
+ */
+#ifdef IEEE_754
+ if(ISNAN(x) )
+ return x + 1.0 ;
+#endif
+ p = x ;
+ if(!R_FINITE(p))
+ return (x < 0.0 ) ? R_DT_0 : R_DT_1;
+ x = p;
+
+ ::Rf_pnorm_both(x, &p, &cp, (lower_tail ? 0 : 1), log_p);
+
+ return(lower_tail ? p : cp);
}
-template <bool NA, typename T>
-inline stats::QNorm<NA,T> qnorm( const Rcpp::VectorBase<REALSXP,NA,T>& x, double mu, double sigma, bool lower = true, bool log = false ) {
- return stats::QNorm<NA,T>( x, mu, sigma, lower, log );
+inline double qnorm_1(double p, double mu /*, double sigma [=1.] */, int lower_tail, int log_p){
+ return ::Rf_qnorm5(p, mu, 1.0, lower_tail, log_p ) ;
}
-
+inline double qnorm_0(double p /*, double mu [=0.], double sigma [=1.] */, int lower_tail, int log_p){
+ return ::Rf_qnorm5(p, 0.0, 1.0, lower_tail, log_p ) ;
}
+} // stats
+} // Rcpp
+
+RCPP_DPQ_0(norm, Rcpp::stats::dnorm_0, Rcpp::stats::pnorm_0, Rcpp::stats::qnorm_0 )
+RCPP_DPQ_1(norm, Rcpp::stats::dnorm_1, Rcpp::stats::pnorm_1, Rcpp::stats::qnorm_1 )
+RCPP_DPQ_2(norm, ::Rf_dnorm4, ::Rf_pnorm5, ::Rf_qnorm5 )
+
#endif
More information about the Rcpp-commits
mailing list