options(error=recover)
set.seed(0xdadada)

##---------------------------------------------------------------
## bag of balls
##---------------------------------------------------------------
n    = 40
p    = 0.1
cols = c("#0000ffd0", "#ff0000d0")

thePlot = function(..., cex=2)
  plot(..., xaxt="n", yaxt="n", bty="n", xlab="", ylab="", pch=16, cex=cex)

x = rep(seq(-1, 1, length.out=n), each=n)
y = rep(seq(-1, 1, length.out=n), times=n)
wh = which((x*x+y*y)<0.6^2)
pdf(file="110629-brixen-deseq-huber-bag.pdf", width=3, height=3)
par(mai=rep(0,4))
thePlot(x=x[wh], y=y[wh], col=ifelse(p<runif(length(wh)), cols[1], cols[2]))
dev.off()


##---------------------------------------------------------------
## Sample from 100 and 5000
##---------------------------------------------------------------
n = c(50, 1000)
for(s in seq(along=n)) {
  nr = c(4, 12)[s]
  nc = n[s]/nr
  nsamp = 4
  x = rep(1:nc, each=nr)
  y = rep(1:nr, times=nc)
  pdf(file=sprintf("110629-brixen-deseq-huber-sample-%d.pdf", n[s]), width=3, height=3)
  par(mai=c(0.1,0,0,0), mfrow=c(nsamp,1))
  for(i in seq_len(nsamp)) 
    thePlot(x, y, col=ifelse(sort(runif(nr*nc)>p), cols[1], cols[2]), ylim=c(0.5,nr+.5), cex=c(3, 0.75)[s])
  dev.off()
}

##---------------------------------------------------------------
## Sample from 100
##---------------------------------------------------------------