# 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

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

step<-1/modenum
curloc<-0
for (i in 1:modenum){
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)
pos<-infopos
}

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