# hierarchischer zufallskritischer Clustertest nach Oldenbürger (1981)
clustest <- function(tab, wdh=100)
{
# (C) LG
# Lizenz: GPL http://www.gnu.org
# zuerst: 13-06-04, zuletzt: 08-06-05, Tübingen
# Permutationstest
# einseitiges Überprüfen ausreichend
# Überschreitungswahrscheinlichkeit vorab festlegen

# Module einbinden für volle Distanzmatrix und Prüfgrösse
  source("optSchnitt.R")	# für volldis()
  source("clustana.R")		# von H. Oldenbürger für DepfuI()

# bei tab wird erwartet: Zeilen = auszuwerten, Spalten = Variablen
  dm <- volldis(dist(tab, method="euclidean"), both=TRUE)
  depfui.o <- as.data.frame(t(DepfuI(dm)))
  colnames(depfui.o) <- c("Depfui.A","Depfui.B")
  DM <- dim(dm)[2]

# Matrix erstellen mit pro Spalte zufällig gezogenen Werten
# (über IDs der empirischen Distanz-Matrix) _ohne_ zurücklegen (Permutation)
  depfui.s <- data.frame(Depfui.A = rep(0,wdh), Depfui.B = rep(0,wdh))
  for(t in 1:wdh)
  {
   perdm <- matrix(data=NA, ncol=DM, nrow=DM)
   perdm <- apply(perdm, 2, function(x) sample(DM,DM))
   dm <- matrix(dm[perdm],ncol=DM,nrow=DM)
   depfui.s[t,] <- as.data.frame(t(DepfuI(dm)))
  }
return(list("Depfui.orig" = depfui.o,
            "Depfui.sim" = depfui.s))
}

# z-Transfo für bequemen Vergleich simuliert - empirisch
clusz <- function(clust)
{
 depfuia.sim.z <- scale(clust[[2]]$Depfui.A)
 depfuib.sim.z <- scale(clust[[2]]$Depfui.B)
 depfuia.orig.z <- (clust[[1]]$Depfui.A - mean(clust[[2]]$Depfui.A)) /
                    sd(clust[[2]]$Depfui.A)
 depfuib.orig.z <- (clust[[1]]$Depfui.B - mean(clust[[2]]$Depfui.B)) /
                    sd(clust[[2]]$Depfui.B)
return(list("DepfuiA.orig.z" = depfuia.orig.z,
            "DepfuiB.orig.z" = depfuib.orig.z,
            "DepfuiA.sim.z" = depfuia.sim.z,
            "DepfuiB.sm.z" = depfuib.sim.z))
}

# Bsp: Testlauf
# mit Vektor
tab <- round(runif(10)*100)
# alternativ: Matrix
# tab <- 1:100
# dim(tab) <- c(10,10)
dmatrix <- volldis(dist(tab), both=TRUE)
DepfuI(dmatrix)
clusres <- clustest(tab,wdh=100)
clusz(clusres)