Visualising hierarchies

hierarchies.R

# Hierarchical uniform distributions
Nsamps=500000
r_min=1
r_max=10
trunc=FALSE

makePlot=function(){
  xlim=c(-5,15)
  op=par(mfrow=c(1,3))
  hist(r,breaks=100,freq=FALSE,xlim=xlim)
  abline(v=c(r_min,r_max),col="red",lwd=2,lty=2)
  hist(r_gen,breaks=100,freq=FALSE,xlim=xlim)
  abline(v=c(r_min,r_max),col="red",lwd=2,lty=2)
  hist(r_rep,breaks=100,freq=FALSE,xlim=xlim)
  abline(v=c(r_min,r_max),col="red",lwd=2,lty=2)
  par(op)
}

pdf("Hierarchies.pdf",width=12, height=6)

r=runif(Nsamps,r_min,r_max)
r_delta=runif(Nsamps,0,(r_max-r_min)/2)

r_gen=runif(Nsamps,r-r_delta,r+r_delta)
r_delta_gen=runif(Nsamps,0,r_delta)
r_rep=runif(Nsamps,r_gen-r_delta_gen,r_gen+r_delta_gen)

makePlot()

r_gen=runif(Nsamps,pmax(r_min,r-r_delta),pmin(r_max,r+r_delta))
r_delta_gen=runif(Nsamps,0,r_delta)
r_rep=runif(Nsamps,pmax(r_min,r_gen-r_delta_gen),pmin(r_max,r_gen+r_delta_gen))

makePlot()

# Hierarchical Normal & Uniform Distributions
r_mu=(r_max+r_min)/2
r_sigma_low=1
r_sigma_up=3

r_sigma_pop=runif(Nsamps,r_sigma_low,r_sigma_up)
r=rnorm(Nsamps,r_mu,r_sigma_pop)

r_sigma_delta_gene=runif(Nsamps,0,(r_sigma_up-r_sigma_low)/2)
r_sigma_gene=runif(Nsamps,pmax(r_sigma_low,r_sigma_pop-r_sigma_delta_gene),pmin(r_sigma_up,r_sigma_pop+r_sigma_delta_gene))
r_gen=rnorm(Nsamps,r,r_sigma_gene)

r_rep=rnorm(Nsamps,r_gen,r_sigma_gene)

makePlot()

dev.off()