diff --git a/DESCRIPTION b/DESCRIPTION index e1c6e54..0903fc6 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,12 +1,10 @@ Package: smcUtils Type: Package -Title: Utility functions for sequential Monte Carlo -Version: 0.2.3 -Author: Jarad Niemi -Maintainer: Jarad Niemi -Date: 2012-12-27 +Title: Utility Functions for Sequential Monte Carlo +Version: 0.2.6 Author: Jarad Niemi Maintainer: Jarad Niemi +Date: 2019-05-02 Description: Provides resampling functions (stratified, residual, multinomial, systematic, and branching), measures of weight uniformity (coefficient of variation, effective sample size, and entropy), and a weight renormalizing function. @@ -15,4 +13,4 @@ URL: https://github.com/jarad/smcUtils LazyLoad: yes Depends: stats Suggests: testthat, Hmisc - +RoxygenNote: 6.1.1 diff --git a/R/branching.resample.R b/R/branching.resample.R index d152bf4..c135f9c 100644 --- a/R/branching.resample.R +++ b/R/branching.resample.R @@ -4,7 +4,7 @@ function(weights,num.samples=length(weights), engine="R") { expected.num.samples = weights*num.samples deterministic.reps = floor(expected.num.samples) - random.reps = rbinom(num.samples,1,expected.num.samples-deterministic.reps) + random.reps = stats::rbinom(num.samples,1,expected.num.samples-deterministic.reps) return(rep2id(deterministic.reps+random.reps)) } diff --git a/R/cov.weights.R b/R/cov.weights.R index ffbb87e..1c134b9 100644 --- a/R/cov.weights.R +++ b/R/cov.weights.R @@ -7,7 +7,7 @@ cov.weights = function(weights, engine="R") switch(engine, { # R implementation - return(var(weights)/mean(weights)^2) + return(stats::var(weights)/mean(weights)^2) # Could compute this as # return(mean((length(weights)*weights-1)^2)) diff --git a/R/inverse.cdf.weights.R b/R/inverse.cdf.weights.R index d0a64c1..a340bf7 100644 --- a/R/inverse.cdf.weights.R +++ b/R/inverse.cdf.weights.R @@ -1,5 +1,5 @@ -inverse.cdf.weights = function(weights, uniforms=runif(length(weights)), engine="R") +inverse.cdf.weights = function(weights, uniforms=stats::runif(length(weights)), engine="R") { check.weights(weights) check.weights(uniforms) diff --git a/R/multinomial.resample.R b/R/multinomial.resample.R index 663afd4..ea9a257 100644 --- a/R/multinomial.resample.R +++ b/R/multinomial.resample.R @@ -14,7 +14,7 @@ multinomial.resample = function(weights, num.samples=length(weights), engine="R" # so to be consistent with C, use inverse.cdf #return(sort(sample(n, num.samples, replace = TRUE, prob = weights))) - return(inverse.cdf.weights(weights,runif(num.samples),engine="R")) + return(inverse.cdf.weights(weights,stats::runif(num.samples),engine="R")) }, { # C implementation diff --git a/R/stratified.resample.R b/R/stratified.resample.R index c5da615..1a9a372 100644 --- a/R/stratified.resample.R +++ b/R/stratified.resample.R @@ -9,7 +9,7 @@ stratified.resample = function(weights, num.samples=length(weights), engine="R") { lbs = seq(0, by=1/num.samples, length=num.samples) ubs = lbs+1/num.samples - u = runif(num.samples, lbs, ubs) + u = stats::runif(num.samples, lbs, ubs) return(inverse.cdf.weights(weights,u,engine="R")) }, { diff --git a/R/systematic.resample.R b/R/systematic.resample.R index aae1af7..5f4a102 100644 --- a/R/systematic.resample.R +++ b/R/systematic.resample.R @@ -9,7 +9,7 @@ systematic.resample = function(weights, num.samples=length(weights), engine="R") switch(engine, { - u = runif(1,0,1/num.samples)+seq(0,by=1/num.samples,length=num.samples) + u = stats::runif(1,0,1/num.samples)+seq(0,by=1/num.samples,length=num.samples) return(inverse.cdf.weights(weights,u,engine="R")) }, { diff --git a/man/resampling.Rd b/man/resampling.Rd index bcb428a..c08519b 100644 --- a/man/resampling.Rd +++ b/man/resampling.Rd @@ -1,80 +1,80 @@ -\name{Unbiased resampling} -\alias{resampling} -\alias{multinomial.resample} -\alias{residual.resample} -\alias{stratified.resample} -\alias{systematic.resample} -\alias{branching.resample} - -\title{Resampling functions} - -\description{A set of resampling functions with unbiased number of replicates.} - -\usage{ -multinomial.resample(weights, num.samples = length(weights), engine="R") -residual.resample( weights, num.samples = length(weights), engine="R", rrf = "stratified") -stratified.resample( weights, num.samples = length(weights), engine="R") -systematic.resample( weights, num.samples = length(weights), engine="R") -branching.resample( weights, num.samples = length(weights), engine="R") -} - -\arguments{ - \item{weights}{a vector of normalized weights} - \item{num.samples}{the number of samples to return (for `branching.resample', - `num.samples' is the expected number of samples as the actual number is random)} - \item{rrf}{for residual resampling, the resampling function to use on the residual} - \item{engine}{run using "R" or "C" code} -} - -\details{ -'multinomial.resample' samples component i with probability `weights[i]', -repeats this sampling `num.samples' times, and returns indices for the sampled -components. - -'residual.resample' deterministically copies `floor(weights)' number of -each component and then performs `rrf' on the remainder. - -`stratified.resample' draws `num.samples' uniform random variables on the -((i-1)/num.samples,i/num.samples) intervals of (0,1). It then uses the -inverse.cdf.weights function to determine which components to sample. - -`systematic.resample' draws 1 uniform random variable on (0,1/num.samples), -builds a sequence of `num.samples' numbers by sequentially adding `1/num.samples', -and then uses `inverse.cdf.weights' to determine which components to -sample. - -`branching.resample' deterministically copies `floor(weights)' number of components -and then draws another component i with probability equal to the residual for -that component. Note: the actual number of components after resampling is random. -} - -\value{Returns a vector of length `num.samples' with indices for sampled components.} - -\references{Douc, R., Cappe, O., Moulines, E. (2005) Comparison of Resampling -Schemes for Particle Filtering. _Proceedings of the 4th International Symposium -on Image and Signal Processing and Analysis_ - -Carpenter, J., Clifford, P., Fearnhead, P. An improved particle filter for -non-linear problems. _IEEE proceedings - Radar, Sonar and Navigation_ *146*, 2-7 -} - -\author{Jarad Niemi} - -%% ~Make other sections like Warning with \section{Warning }{....} ~ - -\seealso{ -\code{\link{resample}},\code{\link{renormalize}}, -} - -\examples{ -ws = renormalize(runif(10)) -multinomial.resample(ws) -residual.resample(ws,rrf="stratified") -stratified.resample(ws,15) -systematic.resample(ws) -branching.resample(ws) -} - -% Add one or more standard keywords, see file 'KEYWORDS' in the -% R documentation directory. -\keyword{ htest } +\name{Unbiased resampling} +\alias{resampling} +\alias{multinomial.resample} +\alias{residual.resample} +\alias{stratified.resample} +\alias{systematic.resample} +\alias{branching.resample} + +\title{Resampling functions} + +\description{A set of resampling functions with unbiased number of replicates.} + +\usage{ +multinomial.resample(weights, num.samples = length(weights), engine="R") +residual.resample(weights, num.samples = length(weights), engine="R", rrf = "stratified") +stratified.resample( weights, num.samples = length(weights), engine="R") +systematic.resample( weights, num.samples = length(weights), engine="R") +branching.resample( weights, num.samples = length(weights), engine="R") +} + +\arguments{ + \item{weights}{a vector of normalized weights} + \item{num.samples}{the number of samples to return (for `branching.resample', + `num.samples' is the expected number of samples as the actual number is random)} + \item{rrf}{for residual resampling, the resampling function to use on the residual} + \item{engine}{run using "R" or "C" code} +} + +\details{ +'multinomial.resample' samples component i with probability `weights[i]', +repeats this sampling `num.samples' times, and returns indices for the sampled +components. + +'residual.resample' deterministically copies `floor(weights)' number of +each component and then performs `rrf' on the remainder. + +`stratified.resample' draws `num.samples' uniform random variables on the +((i-1)/num.samples,i/num.samples) intervals of (0,1). It then uses the +inverse.cdf.weights function to determine which components to sample. + +`systematic.resample' draws 1 uniform random variable on (0,1/num.samples), +builds a sequence of `num.samples' numbers by sequentially adding `1/num.samples', +and then uses `inverse.cdf.weights' to determine which components to +sample. + +`branching.resample' deterministically copies `floor(weights)' number of components +and then draws another component i with probability equal to the residual for +that component. Note: the actual number of components after resampling is random. +} + +\value{Returns a vector of length `num.samples' with indices for sampled components.} + +\references{Douc, R., Cappe, O., Moulines, E. (2005) Comparison of Resampling +Schemes for Particle Filtering. _Proceedings of the 4th International Symposium +on Image and Signal Processing and Analysis_ + +Carpenter, J., Clifford, P., Fearnhead, P. An improved particle filter for +non-linear problems. _IEEE proceedings - Radar, Sonar and Navigation_ *146*, 2-7 +} + +\author{Jarad Niemi} + +%% ~Make other sections like Warning with \section{Warning }{....} ~ + +\seealso{ +\code{\link{resample}},\code{\link{renormalize}}, +} + +\examples{ +ws = renormalize(runif(10)) +multinomial.resample(ws) +residual.resample(ws,rrf="stratified") +stratified.resample(ws,15) +systematic.resample(ws) +branching.resample(ws) +} + +% Add one or more standard keywords, see file 'KEYWORDS' in the +% R documentation directory. +\keyword{ htest } diff --git a/src/init.c b/src/init.c new file mode 100644 index 0000000..88b2b8b --- /dev/null +++ b/src/init.c @@ -0,0 +1,38 @@ +#include // for NULL +#include + +/* FIXME: + Check these declarations against the C/Fortran source code. + */ + +/* .C calls */ +extern void cov2_R(void *, void *, void *); +extern void entropy_R(void *, void *, void *); +extern void ess_R(void *, void *, void *); +extern void inverse_cdf_weights_R(void *, void *, void *, void *, void *); +extern void multinomial_resample_R(void *, void *, void *, void *); +extern void renormalize_R(void *, void *, void *); +extern void rep2id_R(void *, void *, void *); +extern void residual_resample_R(void *, void *, void *, void *, void *); +extern void stratified_resample_R(void *, void *, void *, void *); +extern void systematic_resample_R(void *, void *, void *, void *); + +static const R_CMethodDef CEntries[] = { + {"cov2_R", (DL_FUNC) &cov2_R, 3}, + {"entropy_R", (DL_FUNC) &entropy_R, 3}, + {"ess_R", (DL_FUNC) &ess_R, 3}, + {"inverse_cdf_weights_R", (DL_FUNC) &inverse_cdf_weights_R, 5}, + {"multinomial_resample_R", (DL_FUNC) &multinomial_resample_R, 4}, + {"renormalize_R", (DL_FUNC) &renormalize_R, 3}, + {"rep2id_R", (DL_FUNC) &rep2id_R, 3}, + {"residual_resample_R", (DL_FUNC) &residual_resample_R, 5}, + {"stratified_resample_R", (DL_FUNC) &stratified_resample_R, 4}, + {"systematic_resample_R", (DL_FUNC) &systematic_resample_R, 4}, + {NULL, NULL, 0} +}; + +void R_init_smcUtils(DllInfo *dll) +{ + R_registerRoutines(dll, CEntries, NULL, NULL, NULL); + R_useDynamicSymbols(dll, FALSE); +} \ No newline at end of file