[r-cran-mcmc] 06/11: New upstream version 0.9-5
Andreas Tille
tille at debian.org
Wed Sep 6 20:05:44 UTC 2017
This is an automated email from the git hooks/post-receive script.
tille pushed a commit to branch master
in repository r-cran-mcmc.
commit 6a34f92545013230690a20c0eed5857d1955e7d7
Author: Andreas Tille <tille at debian.org>
Date: Wed Sep 6 21:58:09 2017 +0200
New upstream version 0.9-5
---
.Rinstignore | 3 +-
ChangeLog | 4 +
DESCRIPTION | 15 ++--
MD5 | 77 +++++++++--------
NAMESPACE | 3 +-
R/initseq.R | 2 +-
R/metrop.R | 8 +-
R/olbm.R | 4 +-
R/temper.R | 20 +++--
build/vignette.rds | Bin 294 -> 293 bytes
inst/designDoc/Makefile | 14 +++
inst/{doc => designDoc}/metrop.tex | 0
inst/{doc => designDoc}/temper.tex | 0
inst/doc/Makefile | 42 ---------
inst/doc/bfst.pdf | Bin 245607 -> 229577 bytes
inst/doc/debug.pdf | Bin 102701 -> 97927 bytes
inst/doc/demo.pdf | Bin 400106 -> 383460 bytes
inst/doc/metrop.pdf | Bin 167867 -> 158170 bytes
inst/doc/morph.pdf | Bin 678050 -> 270950 bytes
inst/doc/temper.pdf | Bin 173320 -> 162536 bytes
man/metrop.Rd | 158 +++++++++++++++++++++++++++++-----
man/olbm.Rd | 6 +-
man/temper.Rd | 76 ++++++++++------
src/Makevars | 1 +
src/getListElement.c | 31 -------
src/getScalarInteger.c | 31 -------
src/getScalarLogical.c | 31 -------
src/init.c | 28 ++++++
src/initseq.c | 1 +
src/isAllFinite.c | 31 -------
src/mcmc.h | 21 +++++
src/metrop.c | 172 ++++++++++++++++---------------------
src/myutil.h | 34 +-------
src/olbm.c | 36 +-------
src/temper.c | 34 +-------
tests/accept-batch.R | 22 +++++
tests/accept-batch.Rout.save | 49 +++++++++++
tests/initseq.R | 2 +-
tests/initseq.Rout.save | 10 +--
tests/logitmat.Rout.save | 16 ++--
tests/temp-par-witch.R | 18 ++--
tests/temp-par.Rout.save | 24 +++---
tests/temp-ser-witch.Rout.save | 16 ++--
tests/temp-ser.Rout.save | 22 ++---
44 files changed, 551 insertions(+), 511 deletions(-)
diff --git a/.Rinstignore b/.Rinstignore
index 97c972a..b52aafa 100644
--- a/.Rinstignore
+++ b/.Rinstignore
@@ -1,2 +1 @@
-doc/Makefile
-doc/.*[.]tex
+designDoc/Makefile
diff --git a/ChangeLog b/ChangeLog
index 4f9fb8d..3c0b99e 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -70,3 +70,7 @@
0.9-4 cleaned up some tests
import from stats, as required by R-3.3.0
+
+ 0.9-5 use registration of native routines described in Sections 5.4 and
+ 6.15 of Writing R Extensions
+ add calls to cmpfun inside metrop and temper
diff --git a/DESCRIPTION b/DESCRIPTION
index fed811b..ae6d8b3 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,12 +1,12 @@
Package: mcmc
-Version: 0.9-4
-Date: 2015-07-16
+Version: 0.9-5
+Date: 2017-04-15
Title: Markov Chain Monte Carlo
Author: Charles J. Geyer <charlie at stat.umn.edu> and Leif T. Johnson
<ltjohnson at google.com>
Maintainer: Charles J. Geyer <charlie at stat.umn.edu>
-Depends: R (>= 2.10.0)
-Imports: stats
+Depends: R (>= 3.0.2)
+Imports: stats, compiler
Suggests: xtable, Iso
ByteCompile: TRUE
Description: Simulates continuous distributions of random vectors using
@@ -14,12 +14,13 @@ Description: Simulates continuous distributions of random vectors using
R function that evaluates the log unnormalized density. Algorithms
are random walk Metropolis algorithm (function metrop), simulated
tempering (function temper), and morphometric random walk Metropolis
- (Johnson and Geyer, Annals of Statistics, 2012, function morph.metrop),
+ (Johnson and Geyer, 2012, <https://doi.org/10.1214/12-AOS1048>,
+ function morph.metrop),
which achieves geometric ergodicity by change of variable.
License: MIT + file LICENSE
URL: http://www.stat.umn.edu/geyer/mcmc/,
https://github.com/cjgeyer/mcmc
NeedsCompilation: yes
-Packaged: 2015-07-16 21:03:00 UTC; geyer
+Packaged: 2017-04-15 22:07:33 UTC; geyer
Repository: CRAN
-Date/Publication: 2015-07-17 00:31:01
+Date/Publication: 2017-04-16 10:23:50 UTC
diff --git a/MD5 b/MD5
index 906dd48..c5923d9 100644
--- a/MD5
+++ b/MD5
@@ -1,54 +1,59 @@
-53308ef80c1bbb400bd0c73adb34e795 *ChangeLog
-6a7bd914e2d9ff47bdf231dd2a0d2f6c *DESCRIPTION
+138382003fae18121d80493a77c746b0 *ChangeLog
+7528d4ed73466557515a2a18e83349eb *DESCRIPTION
66e0aaa6082a6fa8bd0f666e544a98b0 *LICENSE
-701a83055a16667ee4d0492493d5e8c8 *NAMESPACE
-2ab15c9110b4e34f9252ed1df9e35762 *R/initseq.R
-45122218ce71e9342c4bf1e8cd07389f *R/metrop.R
+c81060b27aff47c1bd4ce68b2c23ec6d *NAMESPACE
+155b2ace5476c068fb4814b8425a47ce *R/initseq.R
+71e5a5334d09404178307b3b050e162c *R/metrop.R
f403046228ade5ac29a4f8f0bb91aabf *R/morph.R
4b899785af3a35a3a773946c6f13b4ca *R/morph.metrop.R
-05f5205389ba50052d1d33ef2d457aea *R/olbm.R
-e633e27d06c3546b183413cc481d4ea3 *R/temper.R
-95ae0c6c4dc2289fc7686979fa66f7e3 *build/vignette.rds
+32904f94f9bb3c753c588bf823fe5736 *R/olbm.R
+8b7c70bee3c0a14385f00819329efa01 *R/temper.R
+7e147cf7b06b190e5aef3a22e9b4d4f0 *build/vignette.rds
c13cfd5bf58c446f5aa035cb7611c3ed *data/foo.txt.gz
fd9bd39298983c002e96c9e7373200a8 *data/logit.txt.gz
-3f6f0fb621fa345776a25a7d5f778484 *inst/doc/Makefile
+e9a42dd6281a9fafded00cf37ece0474 *inst/designDoc/Makefile
+bb2ba2d13ff2fc81e79f0af3f9ba8ce1 *inst/designDoc/metrop.tex
+177420ab597fd25ca27435d726117a89 *inst/designDoc/temper.tex
80b096a40b8dc0ba718641000bab7114 *inst/doc/bfst.R
756267219c99d98813a8154e4b7b46f4 *inst/doc/bfst.Rnw
-db205118740bb021a53eddb27d2bba70 *inst/doc/bfst.pdf
+8f5d2cb8fdad3a7b73e267de092bb857 *inst/doc/bfst.pdf
85f1c6719e9b9dcf3b56c6eb016935ee *inst/doc/debug.R
c5d145108c9efcff74744153fc68bfe0 *inst/doc/debug.Rnw
-fd7e3066d561b1f4b962085d46bbdaad *inst/doc/debug.pdf
+9bd8b4c0426d0637cab5201dcdd0c6c9 *inst/doc/debug.pdf
63f74f8a150711d3ba8065b174636ecb *inst/doc/demo.R
53007c6969b412aed8ca356ff9347e5a *inst/doc/demo.Rnw
-30eab1b82f92c9c7a359806124abf285 *inst/doc/demo.pdf
-0d0c2ce0dfaa641aab62c181227b1cf5 *inst/doc/metrop.pdf
-bb2ba2d13ff2fc81e79f0af3f9ba8ce1 *inst/doc/metrop.tex
+84be0228cc6646f9694faf98744abf60 *inst/doc/demo.pdf
+eee3eab8bad54767874ff95db49780c7 *inst/doc/metrop.pdf
a52b909d82f36678ee7921fd12d80595 *inst/doc/morph.R
018df3c12cabc9de877013fe309d7447 *inst/doc/morph.Rnw
-40164a8fcdacea29a7022f5b4cb52143 *inst/doc/morph.pdf
-728d89fe79b4a873cc5fc846b16d399c *inst/doc/temper.pdf
-177420ab597fd25ca27435d726117a89 *inst/doc/temper.tex
+b55989c80ddec56647c8803aaaadff6e *inst/doc/morph.pdf
+bddbd4104ab4dd35d879d8966e43172d *inst/doc/temper.pdf
fede69d8e361ef6548637fe54060f709 *man/foo.Rd
f77daf0c3cd0922e649a4a6c03a362fd *man/initseq.Rd
828a4296b11b5041ca9ab38eb74c8b2b *man/logit.Rd
-d77bbb9faee227d2f0f6c86ac88c7c6a *man/metrop.Rd
+3f701bc767ea7a47c7fd16458cc0071a *man/metrop.Rd
7f8ec860609476a448d615d70ff074e2 *man/morph.Rd
656dc31b07400ce3202b0719692b6fba *man/morph.metrop.Rd
-d0908e5b3c8bf592f6a420d761988836 *man/olbm.Rd
-118350de2df78a1d966a7e096e80fdd0 *man/temper.Rd
-c05d85819c8f2d55c67b86d8bc09c599 *src/getListElement.c
-214e108c068ecadb4b1567383b37d125 *src/getScalarInteger.c
-4313c1050389a49f26fc50f51b723d6a *src/getScalarLogical.c
-be626015b2fd7067a06eb3392417cde7 *src/initseq.c
-2e581f0775cc807c353c723113809601 *src/isAllFinite.c
-a27470783b1c322fd8beefdc26ef9231 *src/metrop.c
-8a271789517b6b21c5b65ae6b7086733 *src/myutil.h
-e5ccb1e6d0a7395518b7ef35ad7824ce *src/olbm.c
-abab96dfa4cd64801719cf7d39a59307 *src/temper.c
+7c7a9b03f5b124f4e7b8199c1df65e0b *man/olbm.Rd
+bdfc90719ce178e29799d71292e50c5d *man/temper.Rd
+a84ff979fcebdfd681f7bacc511a9dc3 *src/Makevars
+7af348c402a68bcf2bcedd2c87ce4aa2 *src/getListElement.c
+6fecff8a940cfa809b45ebb1523c0d87 *src/getScalarInteger.c
+e1c4052ef70efd7ece2693cdda8eb27e *src/getScalarLogical.c
+ca127a858f1408883f270c294e48d442 *src/init.c
+ed445c70ff467c056b2243fdc2ade515 *src/initseq.c
+9df5ac56584e97897a7ba7c10e391ec5 *src/isAllFinite.c
+94f2888f854ec400c97c16127e7e5991 *src/mcmc.h
+a51ff07fbc322a8d506033e199b82d74 *src/metrop.c
+be70313e2880a45c4ff804b31a18aa2a *src/myutil.h
+f67bdf45f0c4f15f6d0479720f8a265a *src/olbm.c
+2df8112d1e652db998b434ca5dcc6ec8 *src/temper.c
+772ede8af6bc2b8753e998f5cf39edc3 *tests/accept-batch.R
+38d56fcfbd5127790b8002affff4ca58 *tests/accept-batch.Rout.save
30579f6813cad807b8c301897da45de6 *tests/circle.R
b1d6d5e55c62aaa9f331659161fceba1 *tests/circle.Rout.save
-dfbc7b09d19a44e9415ddead1ed5bcfe *tests/initseq.R
-cba98dc947958991168123a62e8e2bd8 *tests/initseq.Rout.save
+e7781d93db17c14e65dd9afee5af1dad *tests/initseq.R
+086ea1cdd4092988cfb7f8c43f24a596 *tests/initseq.Rout.save
01f734e085fc75dc37615d76381f9568 *tests/isotropic.R
7b1a10898e9abe34fab48f6c43aa51a3 *tests/isotropic.Rout.save
ced46ef84ed611b440a653bae2b25805 *tests/logit.R
@@ -64,7 +69,7 @@ cbf893997832a1b30d6d3694ebdd6937 *tests/logitfun.Rout.save
7a2c271c18d53312291dee40418e12bc *tests/logitlogidx.R
3b4b09c246712b9397deed8c89dd0f86 *tests/logitlogidx.Rout.save
6b721b8a80f8f732bda17972d95c3fa3 *tests/logitmat.R
-72ac94e093a31f23a6c22a14dccb7611 *tests/logitmat.Rout.save
+b92f632b513b4af1e57843e7547aeacb *tests/logitmat.Rout.save
11e13ee697c0cb513728b89580c1a22f *tests/logitnegidx.R
7bc0f8e352a432dfc574f54505e2e624 *tests/logitnegidx.Rout.save
c6f8f05fa12d1f9ad774438181c4cb5d *tests/logitsub.R
@@ -83,13 +88,13 @@ dbc222dd8db7af358381261431ee0f7c *tests/morphtoo.Rout.save
1d84beebb0113a970da2a3f337c57ff3 *tests/saveseed.Rout.save
eb13e0e12345cdaa38d2625f813b5ac9 *tests/saveseedmorph.R
48d749cc35ad957b44a12cf2d141273e *tests/saveseedmorph.Rout.save
-25a3e6688b5bce4708356b4118923bb1 *tests/temp-par-witch.R
+5b83e074fcecbc7154205d007e49ff00 *tests/temp-par-witch.R
a3794f2948a37ac905ae40116cdf9307 *tests/temp-par.R
-bf3bc0d39d6c81855fb2ba92b86cf12a *tests/temp-par.Rout.save
+9250646eaf3dc450d59e1cb6353ce00d *tests/temp-par.Rout.save
e0cdab5f26ba0548eacfaac87368b5d5 *tests/temp-ser-witch.R
-8d0c33d90c3b24b34c4098b5e4d0afc2 *tests/temp-ser-witch.Rout.save
+bebac06312e39a1df0e94f7f0fdf585f *tests/temp-ser-witch.Rout.save
3348850e6103ab402de362dcdf2b717b *tests/temp-ser.R
-fa0b5361231a056559c55945c98c247f *tests/temp-ser.Rout.save
+40b5272b2e05b56db8e4eb4b03703935 *tests/temp-ser.Rout.save
756267219c99d98813a8154e4b7b46f4 *vignettes/bfst.Rnw
9020ac34c43c9a166b7e6529f3dcedeb *vignettes/bfst1.rda
16bf2cf64bb79d5d3c7d69cd596fbebf *vignettes/bfst2.rda
diff --git a/NAMESPACE b/NAMESPACE
index 0cd5e75..a7d32ff 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,5 +1,5 @@
-useDynLib(mcmc)
+useDynLib(mcmc, .registration = TRUE, .fixes = "C_")
export(metrop)
export(morph)
@@ -18,3 +18,4 @@ S3method(temper, tempering)
S3method(temper, "function")
importFrom(stats, runif)
+importFrom(compiler, cmpfun)
diff --git a/R/initseq.R b/R/initseq.R
index a4fea93..7fd103a 100644
--- a/R/initseq.R
+++ b/R/initseq.R
@@ -1,5 +1,5 @@
initseq <- function(x) {
stopifnot(is.numeric(x))
stopifnot(is.finite(x))
- .Call("initseq", x - mean(x), PACKAGE = "mcmc")
+ .Call(C_initseq, x - mean(x))
}
diff --git a/R/metrop.R b/R/metrop.R
index 4898dae..314ee68 100644
--- a/R/metrop.R
+++ b/R/metrop.R
@@ -31,14 +31,18 @@ metrop.function <- function(obj, initial, nbatch, blen = 1,
{
if (! exists(".Random.seed")) runif(1)
saveseed <- .Random.seed
+ obj <- cmpfun(obj)
func1 <- function(state) obj(state, ...)
+ func1 <- cmpfun(func1)
env1 <- environment(fun = func1)
if (missing(outfun)) {
func2 <- NULL
env2 <- NULL
outfun <- NULL
} else if (is.function(outfun)) {
+ outfun <- cmpfun(outfun)
func2 <- function(state) outfun(state, ...)
+ func2 <- cmpfun(func2)
env2 <- environment(fun = func2)
} else {
func2 <- outfun
@@ -46,8 +50,8 @@ metrop.function <- function(obj, initial, nbatch, blen = 1,
}
out.time <- system.time(
- out <- .Call("metrop", func1, initial, nbatch, blen, nspac,
- scale, func2, debug, env1, env2, PACKAGE = "mcmc")
+ out <- .Call(C_metrop, func1, initial, nbatch, blen, nspac,
+ scale, func2, debug, env1, env2)
)
out$initial.seed <- saveseed
out$final.seed <- .Random.seed
diff --git a/R/olbm.R b/R/olbm.R
index 0de5ee9..58e7b50 100644
--- a/R/olbm.R
+++ b/R/olbm.R
@@ -12,14 +12,14 @@ olbm <- function(x, batch.length, demean = TRUE) {
mean <- double(p)
no.calc.mean <- FALSE
}
- out <- .C("olbm",
+ out <- .C(C_olbm,
x=x,
n=as.integer(n),
p=as.integer(p),
batch.length=as.integer(batch.length),
mean=as.double(mean),
var=matrix(as.double(0), p, p),
- no.calc.mean=as.logical(no.calc.mean), PACKAGE = "mcmc")
+ no.calc.mean=as.logical(no.calc.mean))
return(out$var)
}
diff --git a/R/temper.R b/R/temper.R
index 3a0936c..f9b61e4 100644
--- a/R/temper.R
+++ b/R/temper.R
@@ -6,14 +6,18 @@ UseMethod("temper")
temper.tempering <- function(obj, initial, neighbors, nbatch, blen = 1,
nspac = 1, scale = 1, outfun, debug = FALSE, parallel = FALSE, ...)
{
- if (missing(initial)) initial <- obj$final
- if (missing(neighbors)) neighbors <- obj$neighbors
+ # like metrop, ignore initial argument
+ initial <- obj$final
+ # makes no (at least little?) sense to change neighbor structure
+ neighbors <- obj$neighbors
if (missing(nbatch)) nbatch <- obj$nbatch
if (missing(blen)) blen <- obj$blen
if (missing(nspac)) nspac <- obj$nspac
if (missing(scale)) scale <- obj$scale
if (missing(debug)) debug <- obj$debug
- if (missing(parallel)) parallel <- obj$parallel
+ # makes no sense to change from parallel to serial or vice versa
+ # size and shape of state wouldn't even be the same
+ parallel <- obj$parallel
assign(".Random.seed", obj$final.seed, .GlobalEnv)
if (missing(outfun)) {
if (is.null(obj$outfun)) {
@@ -35,14 +39,18 @@ temper.function <- function(obj, initial, neighbors, nbatch, blen = 1,
{
if (! exists(".Random.seed")) runif(1)
saveseed <- .Random.seed
+ obj <- cmpfun(obj)
func1 <- function(state) obj(state, ...)
+ func1 <- cmpfun(func1)
env1 <- environment(fun = func1)
if (missing(outfun)) {
func2 <- NULL
env2 <- NULL
outfun <- NULL
} else if (is.function(outfun)) {
+ outfun <- cmpfun(outfun)
func2 <- function(state) outfun(state, ...)
+ func2 <- cmpfun(func2)
env2 <- environment(fun = func2)
}
@@ -60,10 +68,10 @@ temper.function <- function(obj, initial, neighbors, nbatch, blen = 1,
}
out.time <- system.time(
- out <- .Call("temper", func1, initial, neighbors, nbatch, blen, nspac,
- scale, func2, debug, parallel, env1, env2, PACKAGE = "mcmc")
+ out <- .Call(C_temper, func1, initial, neighbors, nbatch, blen, nspac,
+ scale, func2, debug, parallel, env1, env2)
)
- result <- structure(c(list(lud = obj, initial = initial,
+ result <- structure(c(list(lud = obj,
neighbors = neighbors, nbatch = nbatch, blen = blen, nspac = nspac,
scale = scale, outfun = outfun, debug = debug, parallel = parallel,
initial.seed = saveseed, final.seed = .Random.seed, time = out.time),
diff --git a/build/vignette.rds b/build/vignette.rds
index bea2a1d..2e9315d 100644
Binary files a/build/vignette.rds and b/build/vignette.rds differ
diff --git a/inst/designDoc/Makefile b/inst/designDoc/Makefile
new file mode 100644
index 0000000..2d9b43f
--- /dev/null
+++ b/inst/designDoc/Makefile
@@ -0,0 +1,14 @@
+
+all : metrop.pdf temper.pdf clean
+
+metrop.pdf : metrop.tex
+ pdflatex metrop.tex
+ pdflatex metrop.tex
+
+temper.pdf : temper.tex
+ pdflatex temper.tex
+ pdflatex temper.tex
+
+clean :
+ rm -f *.dvi *.aux *.log
+
diff --git a/inst/doc/metrop.tex b/inst/designDoc/metrop.tex
similarity index 100%
rename from inst/doc/metrop.tex
rename to inst/designDoc/metrop.tex
diff --git a/inst/doc/temper.tex b/inst/designDoc/temper.tex
similarity index 100%
rename from inst/doc/temper.tex
rename to inst/designDoc/temper.tex
diff --git a/inst/doc/Makefile b/inst/doc/Makefile
deleted file mode 100644
index d05375c..0000000
--- a/inst/doc/Makefile
+++ /dev/null
@@ -1,42 +0,0 @@
-
-all : bfst.pdf demo.pdf metrop.pdf temper.pdf debug.pdf clean
-
-demo.tex : demo.Rnw
- $(R_HOME)/bin/R CMD Sweave demo.Rnw
-
-demo.pdf : demo.tex
- pdflatex demo.tex
- pdflatex demo.tex
-
-metrop.pdf : metrop.tex
- pdflatex metrop.tex
- pdflatex metrop.tex
-
-temper.pdf : temper.tex
- pdflatex temper.tex
- pdflatex temper.tex
-
-debug.tex : debug.Rnw
- $(R_HOME)/bin/R CMD Sweave debug.Rnw
-
-debug.pdf : debug.tex
- pdflatex debug.tex
- pdflatex debug.tex
-
-bfst.tex : bfst.Rnw
- $(R_HOME)/bin/R CMD Sweave bfst.Rnw
-
-bfst.pdf : bfst.tex
- pdflatex bfst.tex
- pdflatex bfst.tex
-
-morph.tex : morph.Rnw
- $(R_HOME)/bin/R CMD Sweave morph.Rnw
-
-morph.pdf : morph.tex
- pdflatex morph.tex
- pdflatex morph.tex
-
-clean :
- rm -f *.dvi *.aux *.log demo-fig* morph-*.pdf Rplots.*
-
diff --git a/inst/doc/bfst.pdf b/inst/doc/bfst.pdf
index 0a0db33..4f21b11 100644
Binary files a/inst/doc/bfst.pdf and b/inst/doc/bfst.pdf differ
diff --git a/inst/doc/debug.pdf b/inst/doc/debug.pdf
index 2b4d79a..e903977 100644
Binary files a/inst/doc/debug.pdf and b/inst/doc/debug.pdf differ
diff --git a/inst/doc/demo.pdf b/inst/doc/demo.pdf
index 396d0f7..6fc707c 100644
Binary files a/inst/doc/demo.pdf and b/inst/doc/demo.pdf differ
diff --git a/inst/doc/metrop.pdf b/inst/doc/metrop.pdf
index 0dc647d..acd0ba7 100644
Binary files a/inst/doc/metrop.pdf and b/inst/doc/metrop.pdf differ
diff --git a/inst/doc/morph.pdf b/inst/doc/morph.pdf
index 7c140fd..94340c2 100644
Binary files a/inst/doc/morph.pdf and b/inst/doc/morph.pdf differ
diff --git a/inst/doc/temper.pdf b/inst/doc/temper.pdf
index e76d13d..3de22b6 100644
Binary files a/inst/doc/temper.pdf and b/inst/doc/temper.pdf differ
diff --git a/man/metrop.Rd b/man/metrop.Rd
index bcbe1b6..72925ec 100644
--- a/man/metrop.Rd
+++ b/man/metrop.Rd
@@ -10,20 +10,33 @@
\usage{
metrop(obj, initial, nbatch, blen = 1, nspac = 1, scale = 1, outfun,
debug = FALSE, ...)
+\method{metrop}{function}(obj, initial, nbatch, blen = 1, nspac = 1,
+ scale = 1, outfun, debug = FALSE, ...)
+\method{metrop}{metropolis}(obj, initial, nbatch, blen = 1, nspac = 1,
+ scale = 1, outfun, debug = FALSE, ...)
}
\arguments{
- \item{obj}{an R function that evaluates the log unnormalized probability
+ \item{obj}{Either an \R function or an object of class \code{"metropolis"}
+ from a previous invocation of this function.
+
+ If a function, it evaluates the log unnormalized probability
density of the desired equilibrium distribution of the Markov chain.
- First argument is the state vector of the Markov chain. Other arguments
- arbitrary and taken from the \code{...} arguments of this function.
- Should return \code{- Inf} for points of the state space having
+ Its first argument is the state vector of the Markov chain. Other
+ arguments arbitrary and taken from the \code{...} arguments of this
+ function.
+ It should return \code{-Inf} for points of the state space having
probability zero under the desired equilibrium distribution.
- Alternatively, an object of class \code{"metropolis"} from a
- previous run can be supplied, in which case any missing arguments
+ See also Details and Warning.
+
+ If an object of class \code{"metropolis"}, any missing arguments
(including the log unnormalized density function) are taken from
- this object (up until version 0.7-2 this was incorrect with respect
- to the \code{debug} argument, now it applies to it too).}
- \item{initial}{a real vector, the initial state of the Markov chain.}
+ this object. Also \code{initial} is ignored and the initial state
+ of the Markov chain is the final state from the run recorded in
+ \code{obj}.
+ }
+ \item{initial}{a real vector, the initial state of the Markov chain.
+ Must be feasible, see Details. Ignored if \code{obj} is of
+ class \code{"metropolis"}.}
\item{nbatch}{the number of batches.}
\item{blen}{the length of batches.}
\item{nspac}{the spacing of iterations that contribute to batches.}
@@ -45,10 +58,10 @@ by Tierney (1994), with multivariate normal proposal
producing a Markov chain with equilibrium distribution having a specified
unnormalized density. Distribution must be continuous. Support of the
distribution is the support of the density specified by argument \code{obj}.
-The initial state must satisfy \code{obj(state, ...) > - Inf}.
+The initial state must satisfy \code{obj(state, ...) > -Inf}.
Description of a complete MCMC analysis (Bayesian logistic regression)
-using this function can be found in the vignette \code{demo}
-(\url{../doc/demo.pdf}).
+using this function can be found in the vignette
+\code{vignette("demo", "mcmc")}.
Suppose the function coded by the log unnormalized function (either
\code{obj} or \code{obj$lud}) is actually a log unnormalized density,
@@ -70,6 +83,8 @@ situation.
if \code{outfun} is a function, otherwise the dimension of
\code{state[outfun]} if that makes sense, and the dimension
of \code{state} when \code{outfun} is missing.}
+ \item{accept.batch}{a vector of length \code{nbatch}, the batch means
+ of the acceptances.}
\item{initial}{value of argument \code{initial}.}
\item{final}{final state of Markov chain.}
\item{initial.seed}{value of \code{.Random.seed} before the run.}
@@ -94,12 +109,17 @@ are supplied as additional arguments to \code{metrop}.
If \code{outfun} is a function, then both it and the log unnormalized
density function can be defined without \ldots arguments \emph{if they
have exactly the same arguments list} and that works fine. Otherwise it
-doesn't work. Start the definitions \code{ludfun <- function(state, foo)}
-and \code{outfun <- function(state, bar)} and you get an error about
-unused arguments. Instead start the definitions
-\code{ludfun <- function(state, foo, \ldots)}
-and \code{outfun <- function(state, bar, \ldots)}, supply
-\code{foo} and \code{bar} as additional arguments to \code{metrop},
+doesn't work. Define these functions by
+\preformatted{
+ludfun <- function(state, foo)
+outfun <- function(state, bar)
+}
+and you get an error about unused arguments. Instead define these functions by
+\preformatted{
+ludfun <- function(state, foo, \ldots)
+outfun <- function(state, bar, \ldots)
+}
+and supply \code{foo} and \code{bar} as additional arguments to \code{metrop},
and that works fine.
In short, the log unnormalized density function and \code{outfun} need
@@ -108,17 +128,105 @@ when \ldots is left out and sometimes it doesn't.
Of course, one can avoid this whole issue by always defining the log
unnormalized density function and \code{outfun} to have only one argument
-\code{state} and use global variables (objects in the R global environment) to
+\code{state} and use global variables (objects in the \R global environment) to
specify any other information these functions need to use. That too
-follows the R way. But some people consider that bad programming practice.
+follows the \R way. But some people consider that bad programming practice.
+}
+\section{Philosophy of MCMC}{
+This function follows the philosophy of MCMC explained
+the introductory chapter of the
+\emph{Handbook of Markov Chain Monte Carlo} (Geyer, 2011).
+
+This function automatically does batch means in order to reduce
+the size of output and to enable easy calculation of Monte Carlo standard
+errors (MCSE), which measure error due to the Monte Carlo sampling (not
+error due to statistical sampling --- MCSE gets smaller when you run the
+computer longer, but statistical sampling variability only gets smaller
+when you get a larger data set). All of this is explained in the package
+vignette \code{vignette("demo", "mcmc")} and in Section 1.10 of Geyer (2011).
+
+This function does not apparently
+do \dQuote{burn-in} because this concept does not actually help with MCMC
+(Geyer, 2011, Section 1.11.4) but the re-entrant property of this
+function does allow one to do \dQuote{burn-in} if one wants.
+Assuming \code{ludfun}, \code{start.value}, \code{scale}
+have been already defined
+and are, respectively, an \R function coding the log unnormalized density
+of the target distribution, a valid state of the Markov chain,
+and a useful scale factor,
+\preformatted{
+out <- metrop(ludfun, start.value, nbatch = 1, blen = 1e5, scale = scale)
+out <- metrop(out, nbatch = 100, blen = 1000)
+}
+throws away a run of 100 thousand iterations before doing another run of
+100 thousand iterations that is actually useful for analysis, for example,
+\preformatted{
+apply(out$batch, 2, mean)
+apply(out$batch, 2, sd) / sqrt(out$nbatch)
+}
+give estimates of posterior means and their MCSE assuming the batch length
+(here 1000) was long enough to contain almost all of the signifcant
+autocorrelation (see Geyer, 2011, Section 1.10, for more on MCSE).
+The re-entrant property of this function (the second run starts
+where the first one stops) assures that this is really \dQuote{burn-in}.
+
+The re-entrant property allows one to do very long runs without having to
+do them in one invocation of this function.
+\preformatted{
+out2 <- metrop(out)
+out3 <- metrop(out2)
+batch <- rbind(out$batch, out2$batch, out3$batch)
+}
+produces a result as if the first run had been three times as long.
+}
+\section{Tuning}{
+The \code{scale} argument must be adjusted so that the acceptance rate
+is not too low or too high to get reasonable performance. The rule of
+thumb is that the acceptance rate should be about 25\%.
+But this recommendation (Gelman, et al., 1996) is justified by analysis
+of a toy problem (simulating a spherical multivariate normal distribution)
+for which MCMC is unnecessary. There is no reason to believe this is optimal
+for all problems (if it were optimal, a stronger theorem could be proved).
+Nevertheless, it is clear that at very low acceptance rates the chain makes
+little progress (because in most iterations it does not move) and that at
+very high acceptance rates the chain also makes little progress (because
+unless the log unnormalized density is nearly constant, very high acceptance
+rates can only be achieved by very small values of \code{scale} so the
+steps the chain takes are also very small.
+
+Even in the Gelman, et al. (1996) result, the optimal rate for spherical
+multivariate normal depends on dimension. It is 44\% for \eqn{d = 1}
+and 23\% for \eqn{d = \infty}{d = infinity}.
+Geyer and Thompson (1995) have an example, admittedly for
+simulated tempering (see \code{\link{temper}}) rather than random-walk
+Metropolis, in which no acceptance rate less than 70\% produces an ergodic
+Markov chain. Thus 25\% is merely a rule of thumb. We only know we don't
+want too high or too low. Probably 1\% or 99\% is very inefficient.
}
\references{
+Gelman, A., Roberts, G. O., and Gilks, W. R. (1996)
+Efficient Metropolis jumping rules.
+In \emph{Bayesian Statistics 5: Proceedings of the Fifth Valencia
+ International Meeting}. Edited by J. M. Bernardo,
+ J. O. Berger, A. P. Dawid, and A. F. M. Smith.
+Oxford University Press, Oxford, pp. 599--607.
+
+Geyer, C. J. (2011)
+Introduction to MCMC.
+In \emph{Handbook of Markov Chain Monte Carlo}. Edited by S. P. Brooks,
+A. E. Gelman, G. L. Jones, and X. L. Meng.
+Chapman & Hall/CRC, Boca Raton, FL, pp. 3--48.
+
+Geyer, C. J. and Thompson, E. A. (1995).
+Annealing Markov chain Monte Carlo with applications to ancestral inference.
+\emph{Journal of the American Statistical Association} \bold{90} 909--920.
+
Tierney, L. (1994)
Markov chains for exploring posterior distributions (with discussion).
\emph{Annals of Statistics} \bold{22} 1701--1762.
}
\seealso{
-\code{\link{morph.metrop}}
+\code{\link{morph.metrop}} and \code{\link{temper}}
}
\examples{
h <- function(x) if (all(x >= 0) && sum(x) <= 1) return(1) else return(-Inf)
@@ -127,6 +235,7 @@ out$accept
# acceptance rate too low
out <- metrop(out, scale = 0.1)
out$accept
+t.test(out$accept.batch)$conf.int
# acceptance rate o. k. (about 25 percent)
plot(out$batch[ , 1])
# but run length too short (few excursions from end to end of range)
@@ -134,5 +243,12 @@ out <- metrop(out, nbatch = 1e4)
out$accept
plot(out$batch[ , 1])
hist(out$batch[ , 1])
+acf(out$batch[ , 1], lag.max = 250)
+# looks like batch length of 250 is perhaps OK
+out <- metrop(out, blen = 250, nbatch = 100)
+apply(out$batch, 2, mean) # Monte Carlo estimates of means
+apply(out$batch, 2, sd) / sqrt(out$nbatch) # Monte Carlo standard errors
+t.test(out$accept.batch)$conf.int
+acf(out$batch[ , 1]) # appears that blen is long enough
}
\keyword{misc}
diff --git a/man/olbm.Rd b/man/olbm.Rd
index 49362bf..9de3c32 100644
--- a/man/olbm.Rd
+++ b/man/olbm.Rd
@@ -29,11 +29,11 @@ h <- function(x) if (all(x >= 0) && sum(x) <= 1) return(1) else return(-Inf)
out <- metrop(h, rep(0, 5), 1000)
out <- metrop(out, scale = 0.1)
out <- metrop(out, nbatch = 1e4)
-olbm(out$batch, 150)
+foo <- olbm(out$batch, 150)
# monte carlo estimates (true means are same by symmetry)
-apply(out$batch, 1, mean)
+apply(out$batch, 2, mean)
# monte carlo standard errors (true s. d. are same by symmetry)
-sqrt(diag(olbm(out$batch, 150)))
+sqrt(diag(foo))
# check that batch length is reasonable
acf(out$batch, lag.max = 200)
}
diff --git a/man/temper.Rd b/man/temper.Rd
index d3a3f38..d6de4ff 100644
--- a/man/temper.Rd
+++ b/man/temper.Rd
@@ -4,42 +4,62 @@
\alias{temper.tempering}
\title{Simulated Tempering and Umbrella Sampling}
\description{
- Markov chain Monte Carlo for continuous random vectors using parallel
- or serial simulated tempering, also called umbrella sampling. For
- serial tempering the state of the Markov chain is a pair \eqn{(i, x)},
+ Markov chain Monte Carlo (MCMC) for continuous random vectors using
+ parallel or serial tempering, the latter also called umbrella
+ sampling and simulated tempering.
+ The chain simulates \eqn{k} different distributions on the
+ same state space. In parallel tempering, all the distributions are
+ simulated in each iteration. In serial tempering, only one of the
+ distributions is simulated (a random one). In parallel tempering,
+ the state is a \eqn{k \times p}{k by p} matrix, each row of which
+ is the state for one of the distributions.
+ In serial tempering the state of the Markov chain is a pair \eqn{(i, x)},
where \eqn{i} is an integer between 1 and \eqn{k} and \eqn{x} is a vector
of length \eqn{p}. This pair is represented as a single real vector
- \code{c(i, x)}. For parallel tempering the state of the Markov chain
- is vector of vectors \eqn{(x_1, \ldots, x_k)}{(x[1], \ldots, x[k])},
- where each \code{x} is of length \eqn{p}. This vector of vectors is
- represented as a \eqn{k \times p}{k by p} matrix.
+ \code{c(i, x)}. The variable \eqn{i} says which distribution \eqn{x}
+ is a simulation for.
}
\usage{
temper(obj, initial, neighbors, nbatch, blen = 1, nspac = 1, scale = 1,
outfun, debug = FALSE, parallel = FALSE, \dots)
+\method{temper}{function}(obj, initial, neighbors, nbatch,
+ blen = 1, nspac = 1, scale = 1,
+ outfun, debug = FALSE, parallel = FALSE, \dots)
+\method{temper}{tempering}(obj, initial, neighbors, nbatch,
+ blen = 1, nspac = 1, scale = 1,
+ outfun, debug = FALSE, parallel = FALSE, \dots)
}
\arguments{
- \item{obj}{either an \R function or an object of class \code{"tempering"} from
- a previous run. If a function, should evaluate the log unnormalized
+ \item{obj}{either an \R function or an object of class \code{"tempering"}
+ from a previous run.
+
+ If a function, it should evaluate the log unnormalized
density \eqn{\log h(i, x)}{log h(i, x)} of the desired equilibrium
distribution of the Markov chain for serial tempering (the same function
- is used for both serial and parallel tempering, see details below for
- further explanation). If an object, the log unnormalized density function
+ is used for both serial and parallel tempering, see Details below for
+ further explanation).
+
+ If an object of class \code{"tempering"},
+ the log unnormalized density function
is \code{obj$lud}, and missing arguments of \code{temper} are
obtained from the corresponding elements of \code{obj}.
+
The first argument of the log unnormalized density function is the
- state for simulated tempering \eqn{(i, x)} is supplied as an \R vector
- \code{c(i, x)}; other arguments are arbitrary and taken from
+ is an \R vector \code{c(i, x)}, where \code{i} says which distribution
+ \code{x} is supposed to be a simulation for.
+ Other arguments are arbitrary and taken from
the \code{\dots} arguments of \code{temper}. The log unnormalized density
- function should return \code{-Inf} for points of the state space having
+ function should return \code{-Inf} in regions of the state space having
probability zero.}
\item{initial}{for serial tempering, a real vector \code{c(i, x)} as
described above. For parallel tempering, a real
\eqn{k \times p}{k by p} matrix as described above. In either case,
- the initial state of the Markov chain.}
+ the initial state of the Markov chain.
+ Ignored if \code{obj} has class \code{"tempering"}.}
\item{neighbors}{a logical symmetric matrix of dimension \code{k}
by \code{k}. Elements that are \code{TRUE} indicate jumps
- or swaps attempted by the Markov chain.}
+ or swaps attempted by the Markov chain.
+ Ignored if \code{obj} has class \code{"tempering"}.}
\item{nbatch}{the number of batches.}
\item{blen}{the length of batches.}
\item{nspac}{the spacing of iterations that contribute to batches.}
@@ -64,7 +84,8 @@ temper(obj, initial, neighbors, nbatch, blen = 1, nspac = 1, scale = 1,
indicating the current component are returned.}
\item{debug}{if \code{TRUE} extra output useful for testing.}
\item{parallel}{if \code{TRUE} does parallel tempering, if \code{FALSE} does
- serial tempering.}
+ serial tempering.
+ Ignored if \code{obj} has class \code{"tempering"}.}
\item{...}{additional arguments for \code{obj} or \code{outfun}.}
}
\details{
@@ -159,12 +180,17 @@ are supplied as additional arguments to \code{temper} and that works too.
If \code{outfun} is a function, then both it and the log unnormalized
density function can be defined without \ldots arguments \emph{if they
have exactly the same arguments list} and that works fine. Otherwise it
-doesn't work. Start the definitions \code{ludfun <- function(state, foo)}
-and \code{outfun <- function(state, bar)} and you get an error about
-unused arguments. Instead start the definitions
-\code{ludfun <- function(state, foo, \ldots)}
-and \code{outfun <- function(state, bar, \ldots)}, supply
-\code{foo} and \code{bar} as additional arguments to \code{temper},
+doesn't work. Define these functions by
+\preformatted{
+ludfun <- function(state, foo)
+outfun <- function(state, bar)
+}
+and you get an error about unused arguments. Instead define these functions by
+\preformatted{
+ludfun <- function(state, foo, \ldots)
+outfun <- function(state, bar, \ldots)
+}
+and supply \code{foo} and \code{bar} as additional arguments to \code{metrop},
and that works fine.
In short, the log unnormalized density function and \code{outfun} need
@@ -173,9 +199,9 @@ when \ldots is left out and sometimes it doesn't.
Of course, one can avoid this whole issue by always defining the log
unnormalized density function and \code{outfun} to have only one argument
-\code{state} and use global variables (objects in the R global environment) to
+\code{state} and use global variables (objects in the \R global environment) to
specify any other information these functions need to use. That too
-follows the R way. But some people consider that bad programming practice.
+follows the \R way. But some people consider that bad programming practice.
}
\examples{
d <- 9
diff --git a/src/Makevars b/src/Makevars
new file mode 100644
index 0000000..0bb27b7
--- /dev/null
+++ b/src/Makevars
@@ -0,0 +1 @@
+PKG_CFLAGS=$(C_VISIBILITY)
diff --git a/src/getListElement.c b/src/getListElement.c
index a1ef1a4..ee3953f 100644
--- a/src/getListElement.c
+++ b/src/getListElement.c
@@ -1,35 +1,4 @@
-/*
-*
-* mcmc and MCMC package for R
-* Copyright (c) 2005 Charles J. Geyer
-*
-* All rights reserved.
-*
-* Permission is hereby granted, free of charge, to any person obtaining a copy
-* of this software and associated documentation files (the "Software"), to deal
-* in the Software without restriction, including without limitation the rights
-* to use, copy, modify, merge, publish, distribute, and/or sell copies of the
-* Software, and to permit persons to whom the Software is furnished to do so,
-* provided that the above copyright notice(s) and this permission notice appear
-* in all copies of the Software and that both the above copyright notice(s) and
-* this permission notice appear in supporting documentation.
-*
-* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
-* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
-* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT OF THIRD PARTY RIGHTS.
-* IN NO EVENT SHALL THE COPYRIGHT HOLDER OR HOLDERS INCLUDED IN THIS NOTICE
-* BE LIABLE FOR ANY CLAIM, OR ANY SPECIAL INDIRECT OR CONSEQUENTIAL DAMAGES,
-* OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS,
-* WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
-* ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
-*
-* Except as contained in this notice, the name of a copyright holder shall
-* not be used in advertising or otherwise to promote the sale, use or other
-* dealings in this Software without prior written authorization of the
-* copyright holder.
-*/
-
#include <R.h>
#include <Rinternals.h>
#include "myutil.h"
diff --git a/src/getScalarInteger.c b/src/getScalarInteger.c
index 2195630..309251f 100644
--- a/src/getScalarInteger.c
+++ b/src/getScalarInteger.c
@@ -1,35 +1,4 @@
-/*
-*
-* mcmc and MCMC package for R
-* Copyright (c) 2005 Charles J. Geyer
-*
-* All rights reserved.
-*
-* Permission is hereby granted, free of charge, to any person obtaining a copy
-* of this software and associated documentation files (the "Software"), to deal
-* in the Software without restriction, including without limitation the rights
-* to use, copy, modify, merge, publish, distribute, and/or sell copies of the
-* Software, and to permit persons to whom the Software is furnished to do so,
-* provided that the above copyright notice(s) and this permission notice appear
-* in all copies of the Software and that both the above copyright notice(s) and
-* this permission notice appear in supporting documentation.
-*
-* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
-* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
-* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT OF THIRD PARTY RIGHTS.
-* IN NO EVENT SHALL THE COPYRIGHT HOLDER OR HOLDERS INCLUDED IN THIS NOTICE
-* BE LIABLE FOR ANY CLAIM, OR ANY SPECIAL INDIRECT OR CONSEQUENTIAL DAMAGES,
-* OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS,
-* WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
-* ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
-*
-* Except as contained in this notice, the name of a copyright holder shall
-* not be used in advertising or otherwise to promote the sale, use or other
-* dealings in this Software without prior written authorization of the
-* copyright holder.
-*/
-
#include <R.h>
#include <Rinternals.h>
#include "myutil.h"
diff --git a/src/getScalarLogical.c b/src/getScalarLogical.c
index 79e844d..b75f98f 100644
--- a/src/getScalarLogical.c
+++ b/src/getScalarLogical.c
@@ -1,35 +1,4 @@
-/*
-*
-* mcmc and MCMC package for R
-* Copyright (c) 2005 Charles J. Geyer
-*
-* All rights reserved.
-*
-* Permission is hereby granted, free of charge, to any person obtaining a copy
-* of this software and associated documentation files (the "Software"), to deal
-* in the Software without restriction, including without limitation the rights
-* to use, copy, modify, merge, publish, distribute, and/or sell copies of the
-* Software, and to permit persons to whom the Software is furnished to do so,
-* provided that the above copyright notice(s) and this permission notice appear
-* in all copies of the Software and that both the above copyright notice(s) and
-* this permission notice appear in supporting documentation.
-*
-* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
-* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
-* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT OF THIRD PARTY RIGHTS.
-* IN NO EVENT SHALL THE COPYRIGHT HOLDER OR HOLDERS INCLUDED IN THIS NOTICE
-* BE LIABLE FOR ANY CLAIM, OR ANY SPECIAL INDIRECT OR CONSEQUENTIAL DAMAGES,
-* OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS,
-* WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
-* ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
-*
-* Except as contained in this notice, the name of a copyright holder shall
-* not be used in advertising or otherwise to promote the sale, use or other
-* dealings in this Software without prior written authorization of the
-* copyright holder.
-*/
-
#include <R.h>
#include <Rinternals.h>
#include "myutil.h"
diff --git a/src/init.c b/src/init.c
new file mode 100644
index 0000000..459a574
--- /dev/null
+++ b/src/init.c
@@ -0,0 +1,28 @@
+
+#include <Rdefines.h>
+#include <R_ext/Rdynload.h>
+#include <R_ext/Visibility.h>
+#include "mcmc.h"
+
+static R_NativePrimitiveArgType olbm_types[7] =
+ {REALSXP, INTSXP, INTSXP, INTSXP, REALSXP, REALSXP, LGLSXP};
+
+static R_CMethodDef cMethods[] = {
+ {"olbm", (DL_FUNC) &olbm, 7, olbm_types},
+ {NULL, NULL, 0, NULL}
+};
+
+static R_CallMethodDef callMethods[] = {
+ {"metrop", (DL_FUNC) &metrop, 10},
+ {"temper", (DL_FUNC) &temper, 12},
+ {"initseq", (DL_FUNC) &initseq, 1},
+ {NULL, NULL, 0}
+};
+
+void attribute_visible R_init_mcmc(DllInfo *info)
+{
+ R_registerRoutines(info, cMethods, callMethods, NULL, NULL);
+ R_useDynamicSymbols(info, FALSE);
+ R_forceSymbols(info, TRUE);
+}
+
diff --git a/src/initseq.c b/src/initseq.c
index 7fbad2a..7a840f3 100644
--- a/src/initseq.c
+++ b/src/initseq.c
@@ -2,6 +2,7 @@
#include <R.h>
#include <Rinternals.h>
#include "myutil.h"
+#include "mcmc.h"
SEXP initseq(SEXP x)
{
diff --git a/src/isAllFinite.c b/src/isAllFinite.c
index c129b32..f9c0388 100644
--- a/src/isAllFinite.c
+++ b/src/isAllFinite.c
@@ -1,35 +1,4 @@
-/*
-*
-* mcmc and MCMC package for R
-* Copyright (c) 2005 Charles J. Geyer
-*
-* All rights reserved.
-*
-* Permission is hereby granted, free of charge, to any person obtaining a copy
-* of this software and associated documentation files (the "Software"), to deal
-* in the Software without restriction, including without limitation the rights
-* to use, copy, modify, merge, publish, distribute, and/or sell copies of the
-* Software, and to permit persons to whom the Software is furnished to do so,
-* provided that the above copyright notice(s) and this permission notice appear
-* in all copies of the Software and that both the above copyright notice(s) and
-* this permission notice appear in supporting documentation.
-*
-* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
-* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
-* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT OF THIRD PARTY RIGHTS.
-* IN NO EVENT SHALL THE COPYRIGHT HOLDER OR HOLDERS INCLUDED IN THIS NOTICE
-* BE LIABLE FOR ANY CLAIM, OR ANY SPECIAL INDIRECT OR CONSEQUENTIAL DAMAGES,
-* OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS,
-* WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
-* ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
-*
-* Except as contained in this notice, the name of a copyright holder shall
-* not be used in advertising or otherwise to promote the sale, use or other
-* dealings in this Software without prior written authorization of the
-* copyright holder.
-*/
-
#include <R.h>
#include <Rinternals.h>
#include "myutil.h"
diff --git a/src/mcmc.h b/src/mcmc.h
new file mode 100644
index 0000000..d72c57c
--- /dev/null
+++ b/src/mcmc.h
@@ -0,0 +1,21 @@
+
+#ifndef MCMC_MCMC_H
+#define MCMC_MCMC_H
+
+#include <R.h>
+#include <Rinternals.h>
+
+SEXP metrop(SEXP func1, SEXP initial, SEXP nbatch, SEXP blen, SEXP nspac,
+ SEXP scale, SEXP func2, SEXP debug, SEXP rho1, SEXP rho2);
+
+SEXP temper(SEXP func1, SEXP initial, SEXP neighbors, SEXP nbatch,
+ SEXP blen, SEXP nspac, SEXP scale, SEXP func2, SEXP debug,
+ SEXP parallel, SEXP rho1, SEXP rho2);
+
+SEXP initseq(SEXP x);
+
+void olbm(double *x, int *nin, int *pin, int *lin, double *mean,
+ double *var, int *nocalcin);
+
+#endif /* MCMC_MCMC_H */
+
diff --git a/src/metrop.c b/src/metrop.c
index 848d09a..3df1573 100644
--- a/src/metrop.c
+++ b/src/metrop.c
@@ -1,39 +1,9 @@
-/*
-*
-* mcmc and MCMC package for R
-* Copyright (c) 2005 Charles J. Geyer
-*
-* All rights reserved.
-*
-* Permission is hereby granted, free of charge, to any person obtaining a copy
-* of this software and associated documentation files (the "Software"), to deal
-* in the Software without restriction, including without limitation the rights
-* to use, copy, modify, merge, publish, distribute, and/or sell copies of the
-* Software, and to permit persons to whom the Software is furnished to do so,
-* provided that the above copyright notice(s) and this permission notice appear
-* in all copies of the Software and that both the above copyright notice(s) and
-* this permission notice appear in supporting documentation.
-*
-* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
-* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
-* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT OF THIRD PARTY RIGHTS.
-* IN NO EVENT SHALL THE COPYRIGHT HOLDER OR HOLDERS INCLUDED IN THIS NOTICE
-* BE LIABLE FOR ANY CLAIM, OR ANY SPECIAL INDIRECT OR CONSEQUENTIAL DAMAGES,
-* OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS,
-* WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
-* ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
-*
-* Except as contained in this notice, the name of a copyright holder shall
-* not be used in advertising or otherwise to promote the sale, use or other
-* dealings in this Software without prior written authorization of the
-* copyright holder.
-*/
-
#include <R.h>
#include <Rinternals.h>
#include <Rmath.h>
#include "myutil.h"
+#include "mcmc.h"
static void proposal_setup(SEXP scale, int d);
@@ -52,12 +22,10 @@ SEXP metrop(SEXP func1, SEXP initial, SEXP nbatch, SEXP blen, SEXP nspac,
SEXP state, proposal;
int dim_state, dim_out;
SEXP result, resultnames, acceptance_rate, path,
- save_initial, save_final;
+ save_initial, save_final, acceptance_rate_batches;
double *batch_buffer;
SEXP out_buffer;
- int ibatch, jbatch, ispac;
- int i, k;
double acceptances = 0.0;
double tries = 0.0;
@@ -108,11 +76,11 @@ SEXP metrop(SEXP func1, SEXP initial, SEXP nbatch, SEXP blen, SEXP nspac,
PROTECT(out_buffer = allocVector(REALSXP, dim_out));
if (! int_debug) {
- PROTECT(result = allocVector(VECSXP, 4));
- PROTECT(resultnames = allocVector(STRSXP, 4));
+ PROTECT(result = allocVector(VECSXP, 5));
+ PROTECT(resultnames = allocVector(STRSXP, 5));
} else {
- PROTECT(result = allocVector(VECSXP, 10));
- PROTECT(resultnames = allocVector(STRSXP, 10));
+ PROTECT(result = allocVector(VECSXP, 11));
+ PROTECT(resultnames = allocVector(STRSXP, 11));
}
PROTECT(acceptance_rate = allocVector(REALSXP, 1));
SET_VECTOR_ELT(result, 0, acceptance_rate);
@@ -120,33 +88,39 @@ SEXP metrop(SEXP func1, SEXP initial, SEXP nbatch, SEXP blen, SEXP nspac,
SET_VECTOR_ELT(result, 1, path);
PROTECT(save_initial = duplicate(state));
SET_VECTOR_ELT(result, 2, save_initial);
- UNPROTECT(3);
+ /* cannot set final yet because we haven't got it yet
+ (final value at end of run).
+ See third to last statement of this function. */
+ PROTECT(acceptance_rate_batches = allocVector(REALSXP, int_nbatch));
+ SET_VECTOR_ELT(result, 4, acceptance_rate_batches);
+ UNPROTECT(4);
SET_STRING_ELT(resultnames, 0, mkChar("accept"));
SET_STRING_ELT(resultnames, 1, mkChar("batch"));
SET_STRING_ELT(resultnames, 2, mkChar("initial"));
SET_STRING_ELT(resultnames, 3, mkChar("final"));
+ SET_STRING_ELT(resultnames, 4, mkChar("accept.batch"));
if (int_debug) {
SEXP spath, ppath, gpath, upath, zpath, apath;
int nn = int_nbatch * int_blen * int_nspac;
PROTECT(spath = allocMatrix(REALSXP, dim_state, nn));
- SET_VECTOR_ELT(result, 4, spath);
+ SET_VECTOR_ELT(result, 5, spath);
PROTECT(ppath = allocMatrix(REALSXP, dim_state, nn));
- SET_VECTOR_ELT(result, 5, ppath);
+ SET_VECTOR_ELT(result, 6, ppath);
PROTECT(gpath = allocVector(REALSXP, nn));
- SET_VECTOR_ELT(result, 6, gpath);
+ SET_VECTOR_ELT(result, 7, gpath);
PROTECT(upath = allocVector(REALSXP, nn));
- SET_VECTOR_ELT(result, 7, upath);
+ SET_VECTOR_ELT(result, 8, upath);
PROTECT(zpath = allocMatrix(REALSXP, dim_state, nn));
- SET_VECTOR_ELT(result, 8, zpath);
+ SET_VECTOR_ELT(result, 9, zpath);
PROTECT(apath = allocVector(LGLSXP, nn));
- SET_VECTOR_ELT(result, 9, apath);
+ SET_VECTOR_ELT(result, 10, apath);
UNPROTECT(6);
- SET_STRING_ELT(resultnames, 4, mkChar("current"));
- SET_STRING_ELT(resultnames, 5, mkChar("proposal"));
- SET_STRING_ELT(resultnames, 6, mkChar("log.green"));
- SET_STRING_ELT(resultnames, 7, mkChar("u"));
- SET_STRING_ELT(resultnames, 8, mkChar("z"));
- SET_STRING_ELT(resultnames, 9, mkChar("debug.accept"));
+ SET_STRING_ELT(resultnames, 5, mkChar("current"));
+ SET_STRING_ELT(resultnames, 6, mkChar("proposal"));
+ SET_STRING_ELT(resultnames, 7, mkChar("log.green"));
+ SET_STRING_ELT(resultnames, 8, mkChar("u"));
+ SET_STRING_ELT(resultnames, 9, mkChar("z"));
+ SET_STRING_ELT(resultnames, 10, mkChar("debug.accept"));
}
namesgets(result, resultnames);
UNPROTECT(1);
@@ -157,18 +131,19 @@ SEXP metrop(SEXP func1, SEXP initial, SEXP nbatch, SEXP blen, SEXP nspac,
if (current_log_dens == R_NegInf)
error("log unnormalized density -Inf at initial state");
- for (ibatch = 0, k = 0; ibatch < int_nbatch; ibatch++) {
+ for (int ibatch = 0, k = 0; ibatch < int_nbatch; ibatch++) {
- int j;
+ double acceptances_this_batch = 0.0;
+ double tries_this_batch = 0.0;
- for (i = 0; i < dim_out; i++)
+ for (int i = 0; i < dim_out; i++)
batch_buffer[i] = 0.0;
- for (jbatch = 0; jbatch < int_blen; jbatch++) {
+ for (int jbatch = 0; jbatch < int_blen; jbatch++) {
double proposal_log_dens;
- for (ispac = 0; ispac < int_nspac; ispac++) {
+ for (int ispac = 0; ispac < int_nspac; ispac++) {
int accept;
double u = -1.0; /* impossible return from unif_rand() */
@@ -197,14 +172,13 @@ SEXP metrop(SEXP func1, SEXP initial, SEXP nbatch, SEXP blen, SEXP nspac,
if (int_debug) {
int l = ispac + int_nspac * (jbatch + int_blen * ibatch);
int lbase = l * dim_state;
- SEXP spath = VECTOR_ELT(result, 4);
- SEXP ppath = VECTOR_ELT(result, 5);
- SEXP gpath = VECTOR_ELT(result, 6);
- SEXP upath = VECTOR_ELT(result, 7);
- SEXP zpath = VECTOR_ELT(result, 8);
- SEXP apath = VECTOR_ELT(result, 9);
- int lj;
- for (lj = 0; lj < dim_state; lj++) {
+ SEXP spath = VECTOR_ELT(result, 5);
+ SEXP ppath = VECTOR_ELT(result, 6);
+ SEXP gpath = VECTOR_ELT(result, 7);
+ SEXP upath = VECTOR_ELT(result, 8);
+ SEXP zpath = VECTOR_ELT(result, 9);
+ SEXP apath = VECTOR_ELT(result, 10);
+ for (int lj = 0; lj < dim_state; lj++) {
REAL(spath)[lbase + lj] = REAL(state)[lj];
REAL(ppath)[lbase + lj] = REAL(proposal)[lj];
REAL(zpath)[lbase + lj] = z[lj];
@@ -218,24 +192,28 @@ SEXP metrop(SEXP func1, SEXP initial, SEXP nbatch, SEXP blen, SEXP nspac,
}
if (accept) {
- int jj;
- for (jj = 0; jj < dim_state; jj++)
+ for (int jj = 0; jj < dim_state; jj++)
REAL(state)[jj] = REAL(proposal)[jj];
current_log_dens = proposal_log_dens;
acceptances++;
+ acceptances_this_batch++;
}
tries++;
+ tries_this_batch++;
} /* end of inner loop (one iteration) */
outfun(state, out_buffer);
- for (j = 0; j < dim_out; j++)
+ for (int j = 0; j < dim_out; j++)
batch_buffer[j] += REAL(out_buffer)[j];
} /* end of middle loop (one batch) */
- for (j = 0; j < dim_out; j++, k++)
+ for (int j = 0; j < dim_out; j++, k++)
REAL(path)[k] = batch_buffer[j] / int_blen;
+ REAL(acceptance_rate_batches)[ibatch] =
+ acceptances_this_batch / tries_this_batch;
+
} /* end of outer loop */
PutRNGstate();
@@ -290,9 +268,8 @@ static void proposal_setup(SEXP scale, int d)
SEXP bar;
PROTECT(bar = getAttrib(scale, R_DimSymbol));
if (INTEGER(bar)[0] == d && INTEGER(bar)[1] == d) {
- int i;
scale_factor = (double *) R_alloc(d * d, sizeof(double));
- for (i = 0; i < d * d; i++)
+ for (int i = 0; i < d * d; i++)
scale_factor[i] = REAL(foo)[i];
scale_option = FULL;
} else {
@@ -300,9 +277,8 @@ static void proposal_setup(SEXP scale, int d)
}
UNPROTECT(1);
} else if (LENGTH(foo) == d) {
- int i;
scale_factor = (double *) R_alloc(d, sizeof(double));
- for (i = 0; i < d; i++)
+ for (int i = 0; i < d; i++)
scale_factor[i] = REAL(foo)[i];
scale_option = DIAGONAL;
} else if (LENGTH(foo) == 1) {
@@ -318,7 +294,6 @@ static void proposal_setup(SEXP scale, int d)
static void propose(SEXP state, SEXP proposal, double *z)
{
int d = state_dimension;
- int i, j, k;
if (scale_option == 0)
error("attempt to call propose without setup");
@@ -326,27 +301,27 @@ static void propose(SEXP state, SEXP proposal, double *z)
if (LENGTH(state) != d || LENGTH(proposal) != d)
error("State or proposal length different from initialization\n");
- for (j = 0; j < d; j++)
+ for (int j = 0; j < d; j++)
z[j] = norm_rand();
switch (scale_option) {
case CONSTANT:
- for (j = 0; j < d; j++)
+ for (int j = 0; j < d; j++)
REAL(proposal)[j] = REAL(state)[j]
+ scale_factor[0] * z[j];
break;
case DIAGONAL:
- for (j = 0; j < d; j++)
+ for (int j = 0; j < d; j++)
REAL(proposal)[j] = REAL(state)[j]
+ scale_factor[j] * z[j];
break;
case FULL:
- for (j = 0; j < d; j++)
+ for (int j = 0; j < d; j++)
REAL(proposal)[j] = REAL(state)[j];
- for (i = 0, k = 0; i < d; i++) {
+ for (int i = 0, k = 0; i < d; i++) {
double u = z[i];
- for (j = 0; j < d; j++)
+ for (int j = 0; j < d; j++)
REAL(proposal)[j] += scale_factor[k++] * u;
}
break;
@@ -382,28 +357,31 @@ static int out_setup(SEXP func, SEXP rho, SEXP state)
out_env = rho;
out_dimension = LENGTH(eval(lang2(func, state), rho));
} else if (isLogical(func)) {
- int i;
if (LENGTH(func) != out_state_dimension)
error("is.logical(outfun) & (length(outfun) != length(initial))");
out_option = OUT_INDEX;
out_index = (int *) R_alloc(out_state_dimension, sizeof(int));
- for (i = 0, out_dimension = 0; i < out_state_dimension; i++) {
+ out_dimension = 0;
+ for (int i = 0; i < out_state_dimension; i++) {
out_index[i] = LOGICAL(func)[i];
out_dimension += out_index[i];
}
} else if (isNumeric(func)) {
SEXP foo;
- int foolen, i;
int foopos = 0;
int fooneg = 0;
PROTECT(foo = coerceVector(func, REALSXP));
- foolen = LENGTH(foo);
- for (i = 0; i < foolen; i++) {
+ int foolen = LENGTH(foo);
+ for (int i = 0; i < foolen; i++) {
double foodble = REAL(foo)[i];
+ if (ISNAN(foodble))
+ error("NA or NaN index for outfun");
+ if (! R_FINITE(foodble))
+ error("-Inf or Inf index for outfun");
int fooint = foodble;
- int fooabs = fooint > 0 ? fooint : (- fooint);
+ int fooabs = fooint >= 0 ? fooint : (- fooint);
- if (foodble == 0)
+ if (fooint == 0)
error("is.numeric(outfun) & any(outfun == 0)");
if (foodble != fooint)
error("is.numeric(outfun) & any(outfun != as.integer(outfun))");
@@ -422,32 +400,34 @@ static int out_setup(SEXP func, SEXP rho, SEXP state)
out_option = OUT_INDEX;
out_index = (int *) R_alloc(out_state_dimension, sizeof(int));
if (foopos > 0) {
- for (i = 0; i < out_state_dimension; i++)
+ for (int i = 0; i < out_state_dimension; i++)
out_index[i] = FALSE;
- for (i = 0; i < foolen; i++) {
+ for (int i = 0; i < foolen; i++) {
int fooint = REAL(foo)[i];
out_index[fooint - 1] = TRUE;
}
} else /* (fooneg > 0) */ {
- for (i = 0; i < out_state_dimension; i++)
+ for (int i = 0; i < out_state_dimension; i++)
out_index[i] = TRUE;
- for (i = 0; i < foolen; i++) {
+ for (int i = 0; i < foolen; i++) {
int fooint = REAL(foo)[i];
int fooabs = (- fooint);
out_index[fooabs - 1] = FALSE;
}
}
- for (i = 0, out_dimension = 0; i < out_state_dimension; i++)
+ out_dimension = 0;
+ for (int i = 0; i < out_state_dimension; i++)
out_dimension += out_index[i];
UNPROTECT(1);
+ } else {
+ error("outfun must be NULL, a function, a numeric vector,"
+ " or a logical vector");
}
return out_dimension;
}
static void outfun(SEXP state, SEXP buffer)
{
- int j, k;
-
if (out_option == 0)
error("attempt to call outfun without setup");
@@ -460,11 +440,11 @@ static void outfun(SEXP state, SEXP buffer)
switch (out_option) {
case OUT_IDENTITY:
- for (j = 0; j < out_state_dimension; j++)
+ for (int j = 0; j < out_state_dimension; j++)
REAL(buffer)[j] = REAL(state)[j];
break;
case OUT_INDEX:
- for (j = 0, k = 0; j < out_state_dimension; j++)
+ for (int j = 0, k = 0; j < out_state_dimension; j++)
if (out_index[j])
REAL(buffer)[k++] = REAL(state)[j];
break;
@@ -481,7 +461,7 @@ static void outfun(SEXP state, SEXP buffer)
error("outfun returned vector with non-finite element");
if (LENGTH(foo) != out_dimension)
error("outfun return vector length changed from initial");
- for (k = 0; k < out_dimension; k++)
+ for (int k = 0; k < out_dimension; k++)
REAL(buffer)[k] = REAL(foo)[k];
UNPROTECT(3);
}
diff --git a/src/myutil.h b/src/myutil.h
index ed42cd1..a76d540 100644
--- a/src/myutil.h
+++ b/src/myutil.h
@@ -1,34 +1,6 @@
-/*
-*
-* mcmc and MCMC package for R
-* Copyright (c) 2005 Charles J. Geyer
-*
-* All rights reserved.
-*
-* Permission is hereby granted, free of charge, to any person obtaining a copy
-* of this software and associated documentation files (the "Software"), to deal
-* in the Software without restriction, including without limitation the rights
-* to use, copy, modify, merge, publish, distribute, and/or sell copies of the
-* Software, and to permit persons to whom the Software is furnished to do so,
-* provided that the above copyright notice(s) and this permission notice appear
-* in all copies of the Software and that both the above copyright notice(s) and
-* this permission notice appear in supporting documentation.
-*
-* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
-* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
-* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT OF THIRD PARTY RIGHTS.
-* IN NO EVENT SHALL THE COPYRIGHT HOLDER OR HOLDERS INCLUDED IN THIS NOTICE
-* BE LIABLE FOR ANY CLAIM, OR ANY SPECIAL INDIRECT OR CONSEQUENTIAL DAMAGES,
-* OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS,
-* WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
-* ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
-*
-* Except as contained in this notice, the name of a copyright holder shall
-* not be used in advertising or otherwise to promote the sale, use or other
-* dealings in this Software without prior written authorization of the
-* copyright holder.
-*/
+#ifndef MCMC_MYUTIL_H
+#define MCMC_MYUTIL_H
#include <R.h>
#include <Rinternals.h>
@@ -38,3 +10,5 @@ int getScalarInteger(SEXP foo, char *argname);
int getScalarLogical(SEXP foo, char *argname);
int isAllFinite(SEXP foo);
+#endif /* MCMC_MYUTIL_H */
+
diff --git a/src/olbm.c b/src/olbm.c
index 26adefb..d02c666 100644
--- a/src/olbm.c
+++ b/src/olbm.c
@@ -1,36 +1,6 @@
-/*
-*
-* mcmc and MCMC package for R
-* Copyright (c) 2005 Charles J. Geyer
-*
-* All rights reserved.
-*
-* Permission is hereby granted, free of charge, to any person obtaining a copy
-* of this software and associated documentation files (the "Software"), to deal
-* in the Software without restriction, including without limitation the rights
-* to use, copy, modify, merge, publish, distribute, and/or sell copies of the
-* Software, and to permit persons to whom the Software is furnished to do so,
-* provided that the above copyright notice(s) and this permission notice appear
-* in all copies of the Software and that both the above copyright notice(s) and
-* this permission notice appear in supporting documentation.
-*
-* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
-* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
-* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT OF THIRD PARTY RIGHTS.
-* IN NO EVENT SHALL THE COPYRIGHT HOLDER OR HOLDERS INCLUDED IN THIS NOTICE
-* BE LIABLE FOR ANY CLAIM, OR ANY SPECIAL INDIRECT OR CONSEQUENTIAL DAMAGES,
-* OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS,
-* WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
-* ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
-*
-* Except as contained in this notice, the name of a copyright holder shall
-* not be used in advertising or otherwise to promote the sale, use or other
-* dealings in this Software without prior written authorization of the
-* copyright holder.
-*/
-
#include <R.h>
+#include "mcmc.h"
/* overlapping batch means for vector time series
*
@@ -49,8 +19,8 @@
#define X(I,J) x[(I) + n * (J)]
#define VAR(I,J) var[(I) + p * (J)]
-void olbm(double *x, long *nin, long *pin, long *lin, double *mean,
- double *var, long *nocalcin)
+void olbm(double *x, int *nin, int *pin, int *lin, double *mean,
+ double *var, int *nocalcin)
{
int n = nin[0];
int p = pin[0];
diff --git a/src/temper.c b/src/temper.c
index 78f3e21..3e5ed1a 100644
--- a/src/temper.c
+++ b/src/temper.c
@@ -1,39 +1,9 @@
-/*
-*
-* mcmc and MCMC package for R
-* Copyright (c) 2009 Charles J. Geyer
-*
-* All rights reserved.
-*
-* Permission is hereby granted, free of charge, to any person obtaining a copy
-* of this software and associated documentation files (the "Software"), to deal
-* in the Software without restriction, including without limitation the rights
-* to use, copy, modify, merge, publish, distribute, and/or sell copies of the
-* Software, and to permit persons to whom the Software is furnished to do so,
-* provided that the above copyright notice(s) and this permission notice appear
-* in all copies of the Software and that both the above copyright notice(s) and
-* this permission notice appear in supporting documentation.
-*
-* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
-* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
-* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT OF THIRD PARTY RIGHTS.
-* IN NO EVENT SHALL THE COPYRIGHT HOLDER OR HOLDERS INCLUDED IN THIS NOTICE
-* BE LIABLE FOR ANY CLAIM, OR ANY SPECIAL INDIRECT OR CONSEQUENTIAL DAMAGES,
-* OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS,
-* WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
-* ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
-*
-* Except as contained in this notice, the name of a copyright holder shall
-* not be used in advertising or otherwise to promote the sale, use or other
-* dealings in this Software without prior written authorization of the
-* copyright holder.
-*/
-
#include <R.h>
#include <Rinternals.h>
#include <Rmath.h>
#include "myutil.h"
+#include "mcmc.h"
#ifdef BLEAT
#include <stdio.h>
@@ -66,6 +36,8 @@ SEXP temper(SEXP func1, SEXP initial, SEXP neighbors, SEXP nbatch,
if (nrows(neighbors) != ncols(neighbors))
error("argument \"neighbors\" must have same row and column dimension");
int ncomp = nrows(neighbors);
+ if (ncomp <= 1)
+ error("must have at least two components");
for (int i = 0; i < ncomp; i++)
for (int j = 0; j < ncomp; j++)
if (LOGICAL(neighbors)[i + ncomp * j] != LOGICAL(neighbors)[j + ncomp * i])
diff --git a/tests/accept-batch.R b/tests/accept-batch.R
new file mode 100644
index 0000000..64fe2b9
--- /dev/null
+++ b/tests/accept-batch.R
@@ -0,0 +1,22 @@
+
+# new feature batching acceptance rates
+
+ set.seed(42)
+
+ library(mcmc)
+
+ h <- function(x) if (all(x >= 0) && sum(x) <= 1) return(1) else return(-Inf)
+ out <- metrop(h, rep(0, 5), nbatch = 100, blen = 100, scale = 0.1,
+ debug = TRUE)
+
+ all.equal(out$accept, mean(out$accept.batch))
+
+ foo <- matrix(out$debug.accept, nrow = out$blen)
+ bar <- colMeans(foo)
+ all.equal(out$accept.batch, bar)
+
+ options(digits = 4) # try to keep insanity of computer arithmetic under control
+
+ out$accept
+ t.test(out$accept.batch)$conf.int
+
diff --git a/tests/accept-batch.Rout.save b/tests/accept-batch.Rout.save
new file mode 100644
index 0000000..b454423
--- /dev/null
+++ b/tests/accept-batch.Rout.save
@@ -0,0 +1,49 @@
+
+R version 3.3.3 (2017-03-06) -- "Another Canoe"
+Copyright (C) 2017 The R Foundation for Statistical Computing
+Platform: x86_64-pc-linux-gnu (64-bit)
+
+R is free software and comes with ABSOLUTELY NO WARRANTY.
+You are welcome to redistribute it under certain conditions.
+Type 'license()' or 'licence()' for distribution details.
+
+R is a collaborative project with many contributors.
+Type 'contributors()' for more information and
+'citation()' on how to cite R or R packages in publications.
+
+Type 'demo()' for some demos, 'help()' for on-line help, or
+'help.start()' for an HTML browser interface to help.
+Type 'q()' to quit R.
+
+>
+> # new feature batching acceptance rates
+>
+> set.seed(42)
+>
+> library(mcmc)
+>
+> h <- function(x) if (all(x >= 0) && sum(x) <= 1) return(1) else return(-Inf)
+> out <- metrop(h, rep(0, 5), nbatch = 100, blen = 100, scale = 0.1,
++ debug = TRUE)
+>
+> all.equal(out$accept, mean(out$accept.batch))
+[1] TRUE
+>
+> foo <- matrix(out$debug.accept, nrow = out$blen)
+> bar <- colMeans(foo)
+> all.equal(out$accept.batch, bar)
+[1] TRUE
+>
+> options(digits = 4) # try to keep insanity of computer arithmetic under control
+>
+> out$accept
+[1] 0.2257
+> t.test(out$accept.batch)$conf.int
+[1] 0.2124 0.2390
+attr(,"conf.level")
+[1] 0.95
+>
+>
+> proc.time()
+ user system elapsed
+ 0.168 0.020 0.184
diff --git a/tests/initseq.R b/tests/initseq.R
index 3346028..9816145 100644
--- a/tests/initseq.R
+++ b/tests/initseq.R
@@ -16,7 +16,7 @@
Gamma[k] < 0
Gamma[k] <- 0
- out <- .Call("initseq", x - mean(x))
+ out <- .Call(mcmc:::C_initseq, x - mean(x))
names(out)
all.equal(gamma[1], out$gamma0)
diff --git a/tests/initseq.Rout.save b/tests/initseq.Rout.save
index b7ec81e..a25a91e 100644
--- a/tests/initseq.Rout.save
+++ b/tests/initseq.Rout.save
@@ -1,7 +1,7 @@
-R version 3.2.1 (2015-06-18) -- "World-Famous Astronaut"
-Copyright (C) 2015 The R Foundation for Statistical Computing
-Platform: i686-pc-linux-gnu (32-bit)
+R version 3.3.3 (2017-03-06) -- "Another Canoe"
+Copyright (C) 2017 The R Foundation for Statistical Computing
+Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
@@ -34,7 +34,7 @@ Type 'q()' to quit R.
[1] TRUE
> Gamma[k] <- 0
>
-> out <- .Call("initseq", x - mean(x))
+> out <- .Call(mcmc:::C_initseq, x - mean(x))
> names(out)
[1] "gamma0" "Gamma.pos" "Gamma.dec" "Gamma.con" "var.pos" "var.dec"
[7] "var.con"
@@ -84,4 +84,4 @@ Iso 0.0-17
>
> proc.time()
user system elapsed
- 3.096 0.048 3.138
+ 0.688 0.016 0.695
diff --git a/tests/logitmat.Rout.save b/tests/logitmat.Rout.save
index 95e2336..2b32e0b 100644
--- a/tests/logitmat.Rout.save
+++ b/tests/logitmat.Rout.save
@@ -1,7 +1,6 @@
-R version 2.13.1 (2011-07-08)
-Copyright (C) 2011 The R Foundation for Statistical Computing
-ISBN 3-900051-07-0
+R version 3.3.3 (2017-03-06) -- "Another Canoe"
+Copyright (C) 2017 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
@@ -76,10 +75,10 @@ Type 'q()' to quit R.
> out.metro <- metrop(logl, as.numeric(coefficients(out)), 1e2,
+ scale = sally, debug = TRUE)
> names(out.metro)
- [1] "accept" "batch" "initial" "final" "current"
- [6] "proposal" "log.green" "u" "z" "debug.accept"
-[11] "initial.seed" "final.seed" "time" "lud" "nbatch"
-[16] "blen" "nspac" "scale" "debug"
+ [1] "accept" "batch" "initial" "final" "accept.batch"
+ [6] "current" "proposal" "log.green" "u" "z"
+[11] "debug.accept" "initial.seed" "final.seed" "time" "lud"
+[16] "nbatch" "blen" "nspac" "scale" "debug"
>
> niter <- out.metro$nbatch * out.metro$blen * out.metro$nspac
> niter == nrow(out.metro$current)
@@ -164,3 +163,6 @@ Type 'q()' to quit R.
[1] TRUE
>
>
+> proc.time()
+ user system elapsed
+ 0.340 0.028 0.359
diff --git a/tests/temp-par-witch.R b/tests/temp-par-witch.R
index 8bacc38..cae9aed 100644
--- a/tests/temp-par-witch.R
+++ b/tests/temp-par-witch.R
@@ -1,10 +1,14 @@
+ if ((! exists("DEBUG")) || (! identical(DEBUG, TRUE))) DEBUG <- FALSE
+
library(mcmc)
options(digits=4) # avoid rounding differences
set.seed(42)
+ save.initial.seed <- .Random.seed
+
d <- 3
witch.which <- 1 - (1 / 2)^(1 / d) * (1 / 4)^(seq(0, 5) / d)
witch.which
@@ -37,7 +41,7 @@
thetas <- matrix(0, ncomp, d)
out <- temper(ludfun, initial = thetas, neighbors = neighbors, nbatch = 50,
- blen = 13, nspac = 7, scale = 0.3456789, parallel = TRUE)
+ blen = 13, nspac = 7, scale = 0.3456789, parallel = TRUE, debug = DEBUG)
names(out)
@@ -56,17 +60,17 @@
return(as.numeric(bar))
}
- out <- temper(out, outfun = outfun)
+ out2 <- temper(out, outfun = outfun)
- colMeans(out$batch)
- apply(out$batch, 2, sd) / sqrt(out$nbatch)
+ colMeans(out2$batch)
+ apply(out2$batch, 2, sd) / sqrt(out$nbatch)
### try again
- out <- temper(out, blen = 103, outfun = outfun)
+ out3 <- temper(out2, blen = 103)
- foo <- cbind(colMeans(out$batch),
- apply(out$batch, 2, sd) / sqrt(out$nbatch))
+ foo <- cbind(colMeans(out3$batch),
+ apply(out3$batch, 2, sd) / sqrt(out$nbatch))
colnames(foo) <- c("means", "MCSE")
foo
diff --git a/tests/temp-par.Rout.save b/tests/temp-par.Rout.save
index 51f35e2..2d7d823 100644
--- a/tests/temp-par.Rout.save
+++ b/tests/temp-par.Rout.save
@@ -1,7 +1,6 @@
-R version 2.13.1 (2011-07-08)
-Copyright (C) 2011 The R Foundation for Statistical Computing
-ISBN 3-900051-07-0
+R version 3.3.3 (2017-03-06) -- "Another Canoe"
+Copyright (C) 2017 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
@@ -41,7 +40,7 @@ x1 0.3362 0.4256 0.790 0.429672
x2 0.8475 0.4701 1.803 0.071394 .
x3 1.5143 0.4426 3.422 0.000622 ***
---
-Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
+Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
@@ -108,13 +107,13 @@ Number of Fisher Scoring iterations: 5
+ blen = 10, nspac = 5, scale = 0.56789, parallel = TRUE, debug = TRUE)
>
> names(out)
- [1] "lud" "initial" "neighbors" "nbatch"
- [5] "blen" "nspac" "scale" "outfun"
- [9] "debug" "parallel" "initial.seed" "final.seed"
-[13] "time" "batch" "acceptx" "accepti"
-[17] "initial" "final" "which" "unif.which"
-[21] "state" "log.hastings" "unif.hastings" "proposal"
-[25] "acceptd" "norm" "unif.choose" "coproposal"
+ [1] "lud" "neighbors" "nbatch" "blen"
+ [5] "nspac" "scale" "outfun" "debug"
+ [9] "parallel" "initial.seed" "final.seed" "time"
+[13] "batch" "acceptx" "accepti" "initial"
+[17] "final" "which" "unif.which" "state"
+[21] "log.hastings" "unif.hastings" "proposal" "acceptd"
+[25] "norm" "unif.choose" "coproposal"
>
> ### check decision about within-component or jump/swap
>
@@ -378,3 +377,6 @@ Number of Fisher Scoring iterations: 5
[1] TRUE
>
>
+> proc.time()
+ user system elapsed
+ 1.376 0.008 1.376
diff --git a/tests/temp-ser-witch.Rout.save b/tests/temp-ser-witch.Rout.save
index 581194e..58aa856 100644
--- a/tests/temp-ser-witch.Rout.save
+++ b/tests/temp-ser-witch.Rout.save
@@ -1,7 +1,6 @@
-R version 2.13.1 (2011-07-08)
-Copyright (C) 2011 The R Foundation for Statistical Computing
-ISBN 3-900051-07-0
+R version 3.3.3 (2017-03-06) -- "Another Canoe"
+Copyright (C) 2017 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
@@ -58,10 +57,10 @@ Type 'q()' to quit R.
+ nbatch = 50, blen = 13, nspac = 7, scale = 0.3456789)
>
> names(out)
- [1] "lud" "initial" "neighbors" "nbatch" "blen"
- [6] "nspac" "scale" "outfun" "debug" "parallel"
-[11] "initial.seed" "final.seed" "time" "batch" "acceptx"
-[16] "accepti" "initial" "final" "ibatch"
+ [1] "lud" "neighbors" "nbatch" "blen" "nspac"
+ [6] "scale" "outfun" "debug" "parallel" "initial.seed"
+[11] "final.seed" "time" "batch" "acceptx" "accepti"
+[16] "initial" "final" "ibatch"
>
> out$acceptx
[1] 0.6388889 0.4385246 0.3631714 0.4885246 0.4709677 0.4735516
@@ -147,3 +146,6 @@ Type 'q()' to quit R.
[1] TRUE
>
>
+> proc.time()
+ user system elapsed
+ 2.740 0.028 2.762
diff --git a/tests/temp-ser.Rout.save b/tests/temp-ser.Rout.save
index f7eece0..0c0482b 100644
--- a/tests/temp-ser.Rout.save
+++ b/tests/temp-ser.Rout.save
@@ -1,7 +1,7 @@
-R version 3.2.1 (2015-06-18) -- "World-Famous Astronaut"
-Copyright (C) 2015 The R Foundation for Statistical Computing
-Platform: i686-pc-linux-gnu (32-bit)
+R version 3.3.3 (2017-03-06) -- "Another Canoe"
+Copyright (C) 2017 The R Foundation for Statistical Computing
+Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
@@ -107,13 +107,13 @@ Number of Fisher Scoring iterations: 5
+ log.pseudo.prior = qux)
>
> names(out)
- [1] "lud" "initial" "neighbors" "nbatch"
- [5] "blen" "nspac" "scale" "outfun"
- [9] "debug" "parallel" "initial.seed" "final.seed"
-[13] "time" "batch" "acceptx" "accepti"
-[17] "initial" "final" "ibatch" "which"
-[21] "unif.which" "state" "log.hastings" "unif.hastings"
-[25] "proposal" "acceptd" "norm" "unif.choose"
+ [1] "lud" "neighbors" "nbatch" "blen"
+ [5] "nspac" "scale" "outfun" "debug"
+ [9] "parallel" "initial.seed" "final.seed" "time"
+[13] "batch" "acceptx" "accepti" "initial"
+[17] "final" "ibatch" "which" "unif.which"
+[21] "state" "log.hastings" "unif.hastings" "proposal"
+[25] "acceptd" "norm" "unif.choose"
>
> apply(out$ibatch, 2, mean)
[1] 0.776 0.170 0.000 0.006 0.024 0.010 0.004 0.010
@@ -375,4 +375,4 @@ Number of Fisher Scoring iterations: 5
>
> proc.time()
user system elapsed
- 5.196 0.036 5.232
+ 1.808 0.020 1.820
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/r-cran-mcmc.git
More information about the debian-science-commits
mailing list