[Rcpp-commits] r1938 - in pkg/Rcpp/inst: include/Rcpp/stats unitTests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Aug 6 05:50:45 CEST 2010
Author: dmbates
Date: 2010-08-06 05:50:44 +0200 (Fri, 06 Aug 2010)
New Revision: 1938
Added:
pkg/Rcpp/inst/include/Rcpp/stats/beta.h
Modified:
pkg/Rcpp/inst/include/Rcpp/stats/stats.h
pkg/Rcpp/inst/unitTests/runit.stats.R
Log:
Added d-p-q beta and tests
Added: pkg/Rcpp/inst/include/Rcpp/stats/beta.h
===================================================================
--- pkg/Rcpp/inst/include/Rcpp/stats/beta.h (rev 0)
+++ pkg/Rcpp/inst/include/Rcpp/stats/beta.h 2010-08-06 03:50:44 UTC (rev 1938)
@@ -0,0 +1,112 @@
+// -*- mode: C++; c-indent-level: 4; c-basic-offset: 4; tab-width: 4 -*-
+//
+// beta.h: Rcpp R/C++ interface class library -- beta distribution
+//
+// Copyright (C) 2010 Dirk Eddelbuettel and Romain Francois
+//
+// This file is part of Rcpp.
+//
+// Rcpp is free software: you can redistribute it and/or modify it
+// under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 2 of the License, or
+// (at your option) any later version.
+//
+// Rcpp is distributed in the hope that it will be useful, but
+// WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with Rcpp. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef Rcpp__stats__beta_h
+#define Rcpp__stats__beta_h
+
+namespace Rcpp {
+
+namespace stats {
+
+namespace impl {
+
+ template <bool NA, typename T>
+ class DBeta : public Rcpp::VectorBase< REALSXP, NA, DBeta<NA,T> > {
+ public:
+ typedef typename Rcpp::VectorBase<REALSXP,NA,T> VEC_TYPE;
+
+ DBeta( const VEC_TYPE& vec_, double a_, double b_, bool log_ = false ) :
+ vec(vec_), a(a_), b(b_), log(log_) {}
+
+ inline double operator[]( int i) const { return ::dbeta( vec[i], a, b, log ); }
+
+ inline int size() const { return vec.size(); }
+
+ private:
+ const VEC_TYPE& vec;
+ double a, b;
+ int log;
+
+ };
+
+ template <bool NA, typename T>
+ class PBeta : public Rcpp::VectorBase< REALSXP, NA, PBeta<NA,T> > {
+ public:
+ typedef typename Rcpp::VectorBase<REALSXP,NA,T> VEC_TYPE;
+
+ PBeta( const VEC_TYPE& vec_, double a_, double b_, bool lowertail_ = true, bool log_ = false ) :
+ vec(vec_), a(a_), b(b_), lowertail(lowertail_), log(log_) {}
+
+ inline double operator[]( int i) const { return ::pbeta( vec[i], a, b, lowertail, log ); }
+
+ inline int size() const { return vec.size(); }
+
+ private:
+ const VEC_TYPE& vec;
+ double a, b;
+ int lowertail, log;
+
+ };
+
+ template <bool NA, typename T>
+ class QBeta : public Rcpp::VectorBase< REALSXP, NA, QBeta<NA,T> > {
+ public:
+ typedef typename Rcpp::VectorBase<REALSXP,NA,T> VEC_TYPE;
+
+ QBeta( const VEC_TYPE& vec_, double a_, double b_, bool lowertail_ = true, bool log_ = false ) :
+ vec(vec_), a(a_), b(b_), lowertail(lowertail_), log(log_) {}
+
+ inline double operator[]( int i) const { return ::qbeta( vec[i], a, b, lowertail, log ); }
+
+ inline int size() const { return vec.size(); }
+
+ private:
+ const VEC_TYPE& vec;
+ double a, b;
+ int lowertail, log;
+
+ };
+
+} // impl
+
+template <bool NA, typename T>
+inline impl::DBeta<NA,T> dbeta( const Rcpp::VectorBase<REALSXP,NA,T>& x, double a, double b, bool log = false ) {
+ return impl::DBeta<NA,T>( x, a, b, log );
+}
+
+template <bool NA, typename T>
+inline impl::PBeta<NA,T> pbeta( const Rcpp::VectorBase<REALSXP,NA,T>& x,
+ double a, double b,
+ bool lowertail = true, bool log = false ) {
+ return impl::PBeta<NA,T>( x, a, b, lowertail, log );
+}
+
+template <bool NA, typename T>
+inline impl::QBeta<NA,T> qbeta( const Rcpp::VectorBase<REALSXP,NA,T>& x,
+ double a, double b,
+ bool lowertail = true, bool log = false ) {
+ return impl::QBeta<NA,T>( x, a, b, lowertail, log );
+}
+
+}
+}
+
+#endif
Modified: pkg/Rcpp/inst/include/Rcpp/stats/stats.h
===================================================================
--- pkg/Rcpp/inst/include/Rcpp/stats/stats.h 2010-08-06 02:22:26 UTC (rev 1937)
+++ pkg/Rcpp/inst/include/Rcpp/stats/stats.h 2010-08-06 03:50:44 UTC (rev 1938)
@@ -23,6 +23,7 @@
#define Rcpp__stats__stats_h
#include <Rcpp/stats/binom.h>
+#include <Rcpp/stats/beta.h>
#include <Rcpp/stats/pois.h>
#include <Rcpp/stats/norm.h>
#include <Rcpp/stats/t.h>
Modified: pkg/Rcpp/inst/unitTests/runit.stats.R
===================================================================
--- pkg/Rcpp/inst/unitTests/runit.stats.R 2010-08-06 02:22:26 UTC (rev 1937)
+++ pkg/Rcpp/inst/unitTests/runit.stats.R 2010-08-06 03:50:44 UTC (rev 1938)
@@ -22,17 +22,29 @@
# definition of all the functions at once
f <- list(
- "runit_dbinom" = list(
- signature( x = "integer" ),
- '
+ "runit_dbeta" = list(
+ signature(x = "numeric",
+ a = "numeric", b = "numeric"),
+ '
+ double aa = as<double>(a), bb = as<double>(b) ;
+ NumericVector xx(x) ;
+ return List::create(
+ _["NoLog"] = stats::dbeta( xx, aa, bb),
+ _["Log"] = stats::dbeta( xx, aa, bb, true )
+ ) ;
+ '),
+
+ "runit_dbinom" = list(
+ signature( x = "integer" ),
+ '
IntegerVector xx(x) ;
return List::create(
_["false"] = stats::dbinom( xx, 10, .5),
_["true"] = stats::dbinom( xx, 10, .5, true )
) ;
- '
- )
- , "runit_dpois" = list(
+ '),
+
+ "runit_dpois" = list(
signature( x = "integer" ),
'
IntegerVector xx(x) ;
@@ -62,6 +74,20 @@
) ;
'
),
+ runit_pbeta = list(
+ signature(x = "numeric",
+ a = "numeric", b = "numeric"),
+ '
+ double aa = as<double>(a), bb = as<double>(b) ;
+ NumericVector xx(x) ;
+ return List::create(
+ _["lowerNoLog"] = stats::pbeta( xx, aa, bb),
+ _["lowerLog"] = stats::pbeta( xx, aa, bb, true, true),
+ _["upperNoLog"] = stats::pbeta( xx, aa, bb, false),
+ _["upperLog"] = stats::pbeta( xx, aa, bb, false, true)
+ ) ;
+ '),
+
"runit_pbinom" = list(
signature( x = "numeric", size = "integer", prob = "numeric" ),
'
@@ -306,7 +332,7 @@
vv <- 0:20
checkEquals(fx(vv),
list(lowerNoLog = ppois(vv, 0.5),
- lowerLog = ppois(vv, 0.5, log=TRUE),
+ lowerLog = ppois(vv, 0.5, log=TRUE),
upperNoLog = ppois(vv, 0.5, lower=FALSE),
upperLog = ppois(vv, 0.5, lower=FALSE, log=TRUE)
),
@@ -322,3 +348,34 @@
),
msg = " stats.qpois.prob")
}
+
+test.stats.dbeta <- function() {
+ fx <- .rcpp.stats$runit_dbeta
+ vv <- seq(0, 1, by = 0.1)
+ a <- 0.5; b <- 2.5
+ checkEquals(fx(vv, a, b),
+ list(NoLog = dbeta(vv, a, b),
+ Log = dbeta(vv, a, b, log=TRUE)
+ ),
+ msg = " stats.qbeta")
+}
+
+
+test.stats.pbeta <- function( ) {
+ fx <- .rcpp.stats$runit_pbeta
+ a <- 0.5; b <- 2.5
+ v <- qbeta(seq(0.0, 1.0, by=0.1), a, b)
+ checkEquals(fx(v, a, b),
+ list(lowerNoLog = pbeta(v, a, b),
+ lowerLog = pbeta(v, a, b, log=TRUE),
+ upperNoLog = pbeta(v, a, b, lower=FALSE),
+ upperLog = pbeta(v, a, b, lower=FALSE, log=TRUE)
+ ),
+ msg = " stats.pbeta" )
+ ## Borrowed from R's d-p-q-r-tests.R
+ x <- c(.01, .10, .25, .40, .55, .71, .98)
+ pbval <- c(-0.04605755624088, -0.3182809860569, -0.7503593555585,
+ -1.241555830932, -1.851527837938, -2.76044482378, -8.149862739881)
+ checkEqualsNumeric(fx(x, 0.8, 2)$upperLog, pbval, msg = " stats.pbeta")
+ checkEqualsNumeric(fx(1-x, 2, 0.8)$lowerLog, pbval, msg = " stats.pbeta")
+}
More information about the Rcpp-commits
mailing list