Chapter 13: Manipulation of density estimates

Figure 1

levalg<-function(tree,
plot=T,data=F,crit=NULL,
modelabel=TRUE,ptext=0,leimat=NULL,symbo=NULL,
info=NULL,infolift=0,infopos=0,infochar=NULL,
xmarginleft=0,xmarginright=0,ymargin=0,
xlim=NULL,ylim=NULL,
col="black",col.axis="black",linecol="black",
pch=21,
dimen=NULL,
binlabel=NULL){    

parent<-tree$parent
level<-tree$level
center<-tree$center
if (is.null(center)){
   nodenum<-length(parent)
   center<-matrix(1,dimen,nodenum)
}

mut<-multitree(parent)    # create multitree 
roots<-mut$roots
child<-mut$child
sibling<-mut$sibling 

if (is.null(dimen)){
  d<-dim(center)[1]
}
else{
  d<-dimen
}
crit<-rep(0,d)            # order so that 1st furthest from origo
sibord<-siborder(mut,crit,center)    

mlkm<-moodilkm(parent)
modloc<-mlkm$modloc   
modenum<-mlkm$lkm  

tree$center<-center
modelinks<-siborToModor(tree)   # make links in right order

itemnum<-length(parent)    
verticalPos<-matrix(0,itemnum,1)

step<-1/modenum
curloc<-0
for (i in 1:modenum){
   curmode<-modelinks[i]   
   verticalPos[curmode]<-curloc
   curloc<-curloc+step
} 


for (i in 1:modenum){
   curnode<-modloc[i]
   par<-parent[curnode]
   while (par>0){
      # calculate mid of children of par
      # go to the end of sibling list
        chi<-child[par]
        summa<-verticalPos[chi]
        childNum<-1
        while(sibling[chi]>0){
           chi<-sibling[chi]
           summa<-summa+verticalPos[chi]
           childNum<-childNum+1
        }                            
        verticalPos[par]<-summa/childNum
        par<-parent[par]
   }
}

vPos<-verticalPos
levve<-level
numero<-0
for (i in 1:length(verticalPos)){
  if (binlabel[i]==1){
     numero<-numero+1
     vPos[numero]<-verticalPos[i]
     levve[numero]<-level[i]
  }
}
vPos<-vPos[1:numero]
levve<-levve[1:numero]


if (is.null(ylim)) ylim<-c(0,max(level)+ptext+ymargin)
xlim<-c(min(verticalPos)-xmarginleft,max(verticalPos)+xmarginright)
plot(vPos,levve,xlab="",ylab="",xlim=xlim,ylim=ylim,xaxt="n",yaxt="n",
col=col,col.axis=col.axis,
pch=pch)  

for (i in 1:itemnum){
    if ((parent[i]>0) && (binlabel[i]==1)){
        xchild<-verticalPos[i]
        ychild<-level[i]
        xparent<-verticalPos[parent[i]]
        yparent<-level[parent[i]]
        segments(xparent,yparent,xchild,ychild,col=linecol)
     }
}                

# lets plot info

if (!is.null(info)){


  ifo<-info
  numero<-0
  for (i in 1:length(info)){
    if (binlabel[i]==1){
       numero<-numero+1
       ifo[numero]<-info[i]
    }
  }
  ifo<-ifo[1:numero]


   nodenum<-numero
   infolocx<-matrix(nodenum,1)
   infolocy<-matrix(nodenum,1)

   for (i in 1:nodenum){
     infolocx[i]<-vPos[i] 
     infolocy[i]<-levve[i]+infolift
   }
   digits<-3
   ifo<-format(ifo,digits=digits)
   adj<-NULL
   pos<-infopos
   text(infolocx,infolocy,ifo,pos,adj)       
}

# lets plot character info

if (!is.null(infochar)){
   nodenum<-itemnum
   infolocx<-matrix(nodenum,1)
   infolocy<-matrix(nodenum,1)
   #
   for (i in 1:nodenum){
     infolocx[i]<-verticalPos[i] 
     infolocy[i]<-level[i]+infolift
   }
   pos<-infopos
   text(infolocx,infolocy,infochar,pos)       
}

# lets plot labels for modes

if (modelabel){

xcoor<-verticalPos
ycoor<-level

mlkm<-moodilkm(parent)
modloc<-mlkm$modloc  
modenum<-length(modloc)
modelocx<-matrix(0,modenum,1)
modelocy<-matrix(0,modenum,1)
if (is.null(leimat)){
   if (is.null(symbo)){
       labels<-paste("M",1:modenum,sep="")
   }
   else{
      labels<-paste(symbo,1:modenum,sep="")
   }
}
else{
   labels<-leimat
}
xcor<-matrix(0,modenum,1)                       
for (i in 1:modenum){
    loc<-modloc[i]
    xcor[i]<-xcoor[loc]
}
modloc<-omaord2(modloc,xcor)
for (i in 1:modenum){
    loc<-modloc[i]
    modelocx[i]<-xcoor[loc]
    modelocy[i]<-ycoor[loc]+ptext
}
text(modelocx,modelocy,labels)      
}
}


siborToModor<-function(tree)
{
data<-plotprof(tree,plot=F,data=T,cutlev=NULL,ptext=0,info=NULL,
infolift=0,infopos=0)
vecs<-data$vecs

parent<-tree$parent
mlkm<-moodilkm(parent)
modloc<-mlkm$modloc
moodinum<-mlkm$lkm    #length(modloc)

xcor<-matrix(0,moodinum,1)
for (i in 1:moodinum){
    loc<-modloc[i]
    xcor[i]<-vecs[loc,1]    
}
modloc<-omaord2(modloc,xcor)       

return(modloc)
}

nodenum<-31
parent<-c(0,1,1,2,2,3,3,4,4,5,5,6,6,7,7,8,8,9,9,10,10,11,11,12,12,13,13,14,14,
15,15)
level<- c(4,3,3,2,2,2,2,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)
volume<-matrix(0,nodenum,1)
tree<-list(parent=parent,level=level,volume=volume)

info<-seq(1,31)

binlabel<-matrix(1,31,1)
binlabel[11]<-0
binlabel[22]<-0
binlabel[23]<-0
binlabel[26]<-0
binlabel[30]<-0

# frame 
  levalg(tree,dimen=1,modelabel=F,info=info,xmarginright=0.05,infopos=-0.5,
  binlabel=binlabel)

Figure 2

# frame 1

plot("","",xlim=c(0,4),ylim=c(0,4),xlab="",ylab="",
pch=21,cex=3)
xgrid0<-c(0,2,4)
ygrid0<-c(0,0,0)
xgrid1<-c(0,2,4)
ygrid1<-c(4,4,4)
segments(xgrid0,ygrid0,xgrid1,ygrid1)

xgrid0<-c(0,0)
ygrid0<-c(0,4)
xgrid1<-c(4,4)
ygrid1<-c(0,4)
segments(xgrid0,ygrid0,xgrid1,ygrid1)
text(1,2,"2")
text(3,2,"3")

# frame 2

plot("","",xlim=c(0,4),ylim=c(0,4),xlab="",ylab="",#axes=FALSE,
pch=21,cex=3)
xgrid0<-c(0,2,4)
ygrid0<-c(0,0,0)
xgrid1<-c(0,2,4)
ygrid1<-c(4,4,4)
segments(xgrid0,ygrid0,xgrid1,ygrid1)

xgrid0<-c(0,0)
ygrid0<-c(0,4)
xgrid1<-c(4,4)
ygrid1<-c(0,4)
segments(xgrid0,ygrid0,xgrid1,ygrid1)

xgrid0<-c(1,3)
ygrid0<-c(0,0)
xgrid1<-c(1,3)
ygrid1<-c(4,4)
segments(xgrid0,ygrid0,xgrid1,ygrid1)

text(0.5,2,"4")
text(1.5,2,"5")
text(2.5,2,"6")
text(3.5,2,"7")

# frame 3

plot("","",xlim=c(0,4),ylim=c(0,4),xlab="",ylab="",#axes=FALSE,
pch=21,cex=3)
xgrid0<-c(0,2,4)
ygrid0<-c(0,0,0)
xgrid1<-c(0,2,4)
ygrid1<-c(4,4,4)
segments(xgrid0,ygrid0,xgrid1,ygrid1)

xgrid0<-c(0,0)
ygrid0<-c(0,4)
xgrid1<-c(4,4)
ygrid1<-c(0,4)
segments(xgrid0,ygrid0,xgrid1,ygrid1)

xgrid0<-c(1,3)
ygrid0<-c(0,0)
xgrid1<-c(1,3)
ygrid1<-c(4,4)
segments(xgrid0,ygrid0,xgrid1,ygrid1)

xgrid0<-c(0)
ygrid0<-c(2)
xgrid1<-c(4)
ygrid1<-c(2)
segments(xgrid0,ygrid0,xgrid1,ygrid1)

text(0.5,1,"8")
text(1.5,1,"10")
text(2.5,1,"12")
text(3.5,1,"14")
text(0.5,3,"9")
text(2.5,3,"13")
text(3.5,3,"15")

# frame 4

plot("","",xlim=c(0,4),ylim=c(0,4),xlab="",ylab="",#axes=FALSE,
pch=21,cex=3)
xgrid0<-c(0,2,4)
ygrid0<-c(0,0,0)
xgrid1<-c(0,2,4)
ygrid1<-c(4,4,4)
segments(xgrid0,ygrid0,xgrid1,ygrid1)

xgrid0<-c(0,0)
ygrid0<-c(0,4)
xgrid1<-c(4,4)
ygrid1<-c(0,4)
segments(xgrid0,ygrid0,xgrid1,ygrid1)

xgrid0<-c(1,3)
ygrid0<-c(0,0)
xgrid1<-c(1,3)
ygrid1<-c(4,4)
segments(xgrid0,ygrid0,xgrid1,ygrid1)

xgrid0<-c(0)
ygrid0<-c(2)
xgrid1<-c(4)
ygrid1<-c(2)
segments(xgrid0,ygrid0,xgrid1,ygrid1)

xgrid0<-c(0,0)
ygrid0<-c(1,3)
xgrid1<-c(4,4)
ygrid1<-c(1,3)
segments(xgrid0,ygrid0,xgrid1,ygrid1)

text(0.5,0.5,"16")
text(1.5,0.5,"20")
text(2.5,0.5,"24")
text(3.5,0.5,"28")
text(0.5,1.5,"17")
text(1.5,1.5,"21")
text(2.5,1.5,"25")
text(3.5,1.5,"29")
text(0.5,2.5,"18")
text(0.5,3.5,"19")
text(2.5,3.5,"27")
text(3.5,3.5,"31")

Figure 3

pcf<-sim.data(N=10,type="1d2modal")

lst<-leafsfirst(pcf)

lift<-0.02
numet<-c("A","B","C","D","E","F","G","H","I","J")

split1<-median(pcf$value)
ala<-(pcf$value=median(pcf$value))
split3<-median(pcf$value[yla])


# frame 1
plot(x="",y="",xlim=c(0,length(pcf$value)),ylim=c(0,max(pcf$value)+lift),
     xlab="",ylab="",cex.axis=1.5)
lev0=0
for (i in 1:length(pcf$value)){
    x1<-i-1
    x2<-i
    lev<-pcf$value[i]
    polygon(c(x1,x2,x2,x1),c(lev0,lev0,lev,lev)) #,col="black",lty="blank")
    text(i-1+0.5,pcf$value[i]+lift,numet[i],cex=1.5)
}

# frame 2
plot(x="",y="",xlim=c(-3.8,4.3),ylim=c(5.5,10),
xlab="",ylab="",xaxt="n",yaxt="n")
points(0,10,cex=2,pch=20)
sii<-0.7
text(0-sii,10,0.10,cex=1.5)
segments(0,10,-2,8)
segments(0,10,2,8)

points(-2,8,cex=2,pch=20)
text(-2-sii,8,0.02,cex=1.5)
text(-2+2*sii,8,labels=c("A,B,H,I,J"),cex=1.5)
segments(-2,8,-3,6)
segments(-2,8,-1,6)

points(2,8,cex=2,pch=20)
text(2-sii,8,0.18,cex=1.5)
text(2+2*sii,8,labels=c("C,D,E,F,G"),cex=1.5)
segments(2,8,1,6)
segments(2,8,3,6)

siii<-0.3
points(-3,6,cex=2,pch=20)
text(-3,6-siii,labels=c("A,J"),cex=1.5)
points(-1,6,cex=2,pch=20)
text(-1,6-siii,labels=c("B,H,I"),cex=1.5)
points(1,6,cex=2,pch=20)
text(1,6-siii,labels=c("D,G"),cex=1.5)
points(3,6,cex=2,pch=20)
text(3,6-siii,labels=c("C,E,F"),cex=1.5)

Figure 4

pcf<-sim.data(N=10,type="1d2modal")

lst<-leafsfirst(pcf)

lift<-0.02
numet<-c("1","2","3","4","5","6","7","8","9","10")
or<-order(pcf$value)
jur<-or
for (i in 1:length(or)) jur[or[i]]<-length(or)-i+1

# frame 1
   plot(x="",y="",xlim=c(0,length(pcf$value)),ylim=c(0,max(pcf$value)+lift),
   xlab="",ylab="")
   lev0=0
   for (i in 1:length(pcf$value)){
      x1<-i-1
      x2<-i
      lev<-pcf$value[i]
      polygon(c(x1,x2,x2,x1),c(lev0,lev0,lev,lev)) 
      text(i-1+0.5,pcf$value[i]+lift,numet[jur[i]])
   }

# frame 2
   info<-numet[length(numet):1]
   plottree(lst,info=info,modelabel=FALSE,infopos=-0.6,xmarginright=0.02)

Figure 5

N<-c(16,16)
et<-sim.data(N=N,type="mulmod")

et$index<-NULL

tt<-leafsfirst(et)

# frame 1
pn<-80
dm<-draw.pcf(et,pnum=c(pn,pn))         
contour(dm$x,dm$y,dm$z)

# frame 2
pn<-40
dm<-draw.pcf(et,pnum=c(pn,pn))         
persp(dm$x,dm$y,dm$z,theta=-50,phi=25, 
xlab="",ylab="",zlab="",ticktype="detailed")

# frame 3
plotvolu(tt,ptext=0.002)

# frame 4
  leafsfirst.visu(tt,et,lkmbound=28)

# frame 5
  leafsfirst.visu(tt,et,lkmbound=57)

# frame 6
  leafsfirst.visu(tt,et,lkmbound=59)

Figure 6

N<-c(16,16)
pcf<-sim.data(N=N,type="fox")
st<-leafsfirst(pcf,propor=0.1)

# frame 1
  leafsfirst.visu(st,pcf,lkmbound=38,propor=0.1)

# frame 2
  leafsfirst.visu(st,pcf,lkmbound=65,propor=0.1)

# frame 3
  leafsfirst.visu(st,pcf,lkmbound=70,propor=0.1)

Figure 7

plotnet<-function(){

plot("","",xlim=c(0,4),ylim=c(0,4),xlab="",ylab="",#axes=FALSE,
pch=21,cex=3)
xgrid0<-c(0,2,4)
ygrid0<-c(0,0,0)
xgrid1<-c(0,2,4)
ygrid1<-c(4,4,4)
segments(xgrid0,ygrid0,xgrid1,ygrid1)

xgrid0<-c(0,0)
ygrid0<-c(0,4)
xgrid1<-c(4,4)
ygrid1<-c(0,4)
segments(xgrid0,ygrid0,xgrid1,ygrid1)

xgrid0<-c(1,3)
ygrid0<-c(0,0)
xgrid1<-c(1,3)
ygrid1<-c(4,4)
segments(xgrid0,ygrid0,xgrid1,ygrid1)

xgrid0<-c(0)
ygrid0<-c(2)
xgrid1<-c(4)
ygrid1<-c(2)
segments(xgrid0,ygrid0,xgrid1,ygrid1)

xgrid0<-c(0,0)
ygrid0<-c(1,3)
xgrid1<-c(4,4)
ygrid1<-c(1,3)
segments(xgrid0,ygrid0,xgrid1,ygrid1)

text(0.5,0.5,"16")
text(1.5,0.5,"20")
text(2.5,0.5,"24")
text(3.5,0.5,"28")
text(0.5,1.5,"17")
text(1.5,1.5,"21")
text(2.5,1.5,"25")
text(3.5,1.5,"29")
text(0.5,2.5,"18")
text(0.5,3.5,"19")
text(2.5,3.5,"27")
text(3.5,3.5,"31")
}

# frame 1

plotnet()

for (i in 0:3){
   x0<-i+0.1
   y0<-0.1
   x1<-i+0.1
   y1<-1.9
   segments(x0,y0,x1,y1)
}
for (i in 0:3){
   x0<-i+0.9
   y0<-0.1
   x1<-i+0.9
   y1<-1.9
   segments(x0,y0,x1,y1)
}
for (i in 0:3){
   x0<-i+0.1
   y0<-0.1
   x1<-i+0.9
   y1<-0.1
   segments(x0,y0,x1,y1)
}
for (i in 0:3){
   x0<-i+0.1
   y0<-1.9
   x1<-i+0.9
   y1<-1.9
   segments(x0,y0,x1,y1)
}

x<-c(0.1,0.1,0.9,0.9,0.1)
y<-c(2.1,3.9,3.9,2.1,2.1)
lines(x,y)

xs<-c(2.1,2.1,2.9,2.9,2.1)
ys<-c(3.1,3.9,3.9,3.1,3.1)
lines(xs,ys)

xu<-xs+1
y<-c(3.1,3.9,3.9,3.1,3.1)
lines(xu,y)

# frame 2

plotnet()

x<-c(0.1,0.1,0.9,0.9,0.1)
y<-c(0.1,3.9,3.9,0.1,0.1)
lines(x,y)

x<-c(1.1,1.1,1.9,1.9,1.1)
y<-c(0.1,1.9,1.9,0.1,0.1)
lines(x,y)

lines(x+1,y)

lines(x+2,y)

xs<-c(2.1,2.1,2.9,2.9,2.1)
ys<-c(3.1,3.9,3.9,3.1,3.1)
lines(xs,ys)

lines(xs+1,ys)

# frame 3

plotnet()

x<-c(0.1,0.1,0.9,0.9,1.9,1.9,0.1)
y<-c(0.1,3.9,3.9,1.9,1.9,0.1,0.1)
lines(x,y)

x<-c(2.1,2.1,3.9,3.9,2.1)
y<-c(0.1,1.9,1.9,0.1,0.1)
lines(x,y)

x<-c(2.1,2.1,3.9,3.9,2.1)
y<-c(3.1,3.9,3.9,3.1,3.1)
lines(x,y)

# frame 4

plotnet()

x<-c(0.1,0.1,0.9,0.9,3.9,3.9,0.1)
y<-c(0.1,3.9,3.9,1.9,1.9,0.1,0.1)
lines(x,y)

x<-c(2.1,2.1,3.9,3.9,2.1)
y<-c(3.1,3.9,3.9,3.1,3.1)
lines(x,y)