Nothing
#### --- Choose cutoff such that i.i.d. re-sample gives MC(0)
library(VLMC)
data(OZrain)
.proctime00 <- proc.time()
rain4c <- cut(OZrain, c(-.1, 0.5, 18.5, 50.1, 1000))
table(rain4c <- as.integer(rain4c) - 1L)
(v4c <- vlmc(rain4c))
##' Finding cutoff for vlmc() such that the order of the fitted Markov chain is practically zero
fff <- function(x, show=getOption("verbose")) {
f <- (s <- vlmc(sr4, cut = x)$size)["ord.MC"]
if(show)
cat(formatC(x,wid=15)," : ", formatC(f,wid=2),
" (",formatC(s["total"],wid=5),")\n", sep="")
f - 0.01
}
N <- 64# had 200
ur.nms <- c("root", "f.root", "iter", "estim.prec")
r <- matrix(NA, length(ur.nms), N,
dimnames = list(ur.nms, NULL))
RNGversion("3.5.0")# + warning .. FIXME once we depend on R >= 3.6.0
set.seed(6352)
for(i in 1:N) {
sr4 <- sample(rain4c)# random permutation -- should be iid!
ur <- uniroot(fff, c(1, 20), tol = 1e-2)
r[,i] <- rr <- unlist(ur[ur.nms])
}
rownames(r) <- names(rr)
r <- t(r)
summary( r[,"root"]) # astonishingly wide; Q1--Q3 = 5.41--6.64
sqrt(var(r[,"root"]))/ sqrt(N) # 0.14
if(!dev.interactive(orNone=TRUE)) pdf("iid-cutoff.pdf")
hist(r[,"root"], xlab = "vlmc - cutoff",
main = paste("min. cutoffs for vlmc(<random perm>) for independence,\n",
"n =", length(rain4c),"; k =", v4c$alpha.len))
## Last Line:
cat('Time elapsed: ', proc.time() - .proctime00,'\n')
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.