-
Notifications
You must be signed in to change notification settings - Fork 22
/
burnin.R
executable file
·54 lines (50 loc) · 1.96 KB
/
burnin.R
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
###########################################################################
# burnin #
# #
# The purpose of the burnin function is to estimate the duration of #
# burn-in in iterations for one or more MCMC chains. #
###########################################################################
burnin <- function(x, method="BMK")
{
if(missing(x)) stop("The x argument is required.")
if(is.vector(x)) x <- matrix(x, length(x), 1)
n <- nrow(x)
burn <- rep(0,ncol(x))
if(method == "BMK") {
if(n %% 10 == 0) x2 <- x
if(n %% 10 != 0) x2 <- x[1:(10*trunc(n/10)),]
HD <- BMK.Diagnostic(x2, 10)
Ind <- 1 * (HD > 0.5)
burn <- n
batch.list <- seq(from=1, to=nrow(x2), by=floor(nrow(x2)/10))
for (i in 1:9) {
if(sum(Ind[,i:9]) == 0) {
burn <- batch.list[i] - 1
break
}
}
}
else {
for (i in 1:ncol(x)) {
iter <- 1
stationary <- 0
jump <- round(n/10)
while(stationary == 0) {
if(method == "KS") {
p <- KS.Diagnostic(x[iter:n])
if(p <= 0.05) stationary <- 1}
else { #method == Geweke
z <- try(Geweke.Diagnostic(x[iter:n]),
silent=TRUE)
if(inherits(z, "try-error")) z <- 3
if(abs(z < 2)) stationary <- 1}
if(stationary == 0) iter <- iter + jump
if(iter >= n) stationary <- 1}
if(iter > 1) iter <- iter - 1
if(iter > n) iter <- n
burn[i] <- iter
}
}
return(burn)
}
#End