File: temper.R

package info (click to toggle)
r-cran-mcmc 0.9-7-1
  • links: PTS, VCS
  • area: main
  • in suites: bookworm, bullseye
  • size: 2,352 kB
  • sloc: ansic: 1,298; makefile: 14; sh: 8
file content (77 lines) | stat: -rw-r--r-- 2,829 bytes parent folder | download | duplicates (2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77

temper <- function(obj, initial, neighbors, nbatch, blen = 1,
    nspac = 1, scale = 1, outfun, debug = FALSE, parallel = FALSE, ...)
UseMethod("temper")

temper.tempering <- function(obj, initial, neighbors, nbatch, blen = 1,
    nspac = 1, scale = 1, outfun, debug = FALSE, parallel = FALSE, ...)
{
    # 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
    # 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)) {
            temper.function(obj$lud, initial, neighbors, nbatch, blen,
                nspac, scale, debug = debug, parallel = parallel, ...)
        } else {
            temper.function(obj$lud, initial, neighbors, nbatch, blen,
                nspac, scale, obj$outfun, debug = debug, parallel = parallel,
                ...)
        }
    } else {
        temper.function(obj$lud, initial, neighbors, nbatch, blen,
            nspac, scale, outfun, debug = debug, parallel = parallel, ...)
    }
}

temper.function <- function(obj, initial, neighbors, nbatch, blen = 1,
    nspac = 1, scale = 1, outfun, debug = FALSE, parallel = FALSE, ...)
{
    if (! exists(".Random.seed")) runif(1)
    saveseed <- .Random.seed
    func1 <- function(state) obj(state, ...)
    env1 <- environment(fun = func1)
    if (missing(outfun)) {
        func2 <- NULL
        env2 <- NULL
        outfun <- NULL
    } else if (is.function(outfun)) {
        func2 <- function(state) outfun(state, ...)
        env2 <- environment(fun = func2)
    }

    stopifnot(is.numeric(initial))
    storage.mode(initial) <- "double"

    if (is.list(scale)) {
        for (i in 1:length(scale)) {
            stopifnot(is.numeric(scale[[i]]))
            storage.mode(scale[[i]]) <- "double"
        }
    } else {
        stopifnot(is.numeric(scale))
        storage.mode(scale) <- "double"
    }

    out.time <- system.time(
    out <- .Call(C_temper, func1, initial, neighbors, nbatch, blen, nspac,
        scale, func2, debug, parallel, env1, env2)
    )
    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),
        out), class = c("mcmc", "tempering"))
    return(result)
}