R Under development (unstable) (2024-08-17 r87027 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2024 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

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.

> #### --- 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)

   0    1    2    3 
2237  713  356  347 
> (v4c <- vlmc(rain4c))

'vlmc' a Variable Length Markov Chain;
	 alphabet '0123', |alphabet| = 4, n = 3653.
Call: vlmc(dts = rain4c);  default cutoff =3.907
 -> extensions (= $size ) :
   ord.MC   context nr.leaves     total 
        9        39        11       322 
AIC =  7485 
> 
> ##' 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
Warning message:
In RNGkind("Mersenne-Twister", "Inversion", "Rounding") :
  non-uniform 'Rounding' sampler used
> 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
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  4.412   5.439   5.996   6.232   6.823  10.494 
> sqrt(var(r[,"root"]))/ sqrt(N) # 0.14
[1] 0.1445622
> 
> 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')
Time elapsed:  10.44 1.85 12.32 NA NA 
> 
> proc.time()
   user  system elapsed 
  10.73    1.95   12.64