R Code for Reproducing the Figures of the LST article.
First issue the command
library(denpro)
library(delt) # for drawing histograms
- 2 Dimensional Data:
Mixture of 3 Gaussians (simulated)
- 3 Dimensional Data:
Mixture of 4 Gaussians (simulated)
- 4 Dimensional Data:
Mixture of 5 Gaussians (simulated)
- Data
- Histogram Estimate
- Kernel Estimate
n<-200
d<-2
moodi<-3
M<-matrix(0,moodi,d)
dist<-4
M[1,]<-dist*c(0,0)
M[2,]<-dist*c(0,1)
M[3,]<-dist*c(1,0)
sig<-matrix(1,moodi,d)
p0<-1/moodi
p<-p0*rep(1,moodi)
seed<-2
dendat<-simmix(n,M,sig,p,seed)
dh8<-drawhist(dendat,binlkm=8,plkm=80)
dh13<-drawhist(dendat,binlkm=13,plkm=80)
ph8<-profhist(dendat,binlkm=8)
ph13<-profhist(dendat,binlkm=13)
postscript(file="~/artlst/histopersp.ps") #,horizontal=FALSE)
par(mfcol=c(2,2))
parold<-par(no.readonly=FALSE)$plt
par(plt=c(0.08,0.95,0.11,0.91))
#par(plt=c(0.05,0.95,0.07,0.93))
persp(dh8$x,dh8$y,dh8$z,col="white",theta=130,phi=20,
xlab="x",ylab="y",zlab="",ticktype="detailed") #,shade=0.5)
text(-0.3,-0.5,"a)",cex=1.5)
#title(sub="a) 8 bins")
persp(dh13$x,dh13$y,dh13$z,col="white",theta=130,phi=20,
xlab="x",ylab="y",zlab="",ticktype="detailed")
text(-0.3,-0.5,"b)",cex=1.5)
#title(sub="b) 13 bins")
par(plt=parold)
#par(plt=c( 0.1172537, 0.9399432, 0.1457143, 0.8828571))
plotvolu(ph8,ptext=0.0025,symbo="L")
plotvolu(ph13,ptext=0.004)
dev.off()
paletti<-c("red","blue","green",
"orange","navy","darkgreen",
"orchid","aquamarine","turquoise",
"pink","violet","magenta","chocolate","cyan",
colors()[50:100])
postscript(file="~/artlst/histobary.ps") #,horizontal=FALSE)
par(mfrow=c(2,2))
col<-colobary(ph8$parent,paletti)
plotbary(ph8,ptext=0.0023,coordi=1,col=col,nodesymbo=19,symbo="L",collines=col)
title(sub="a.1) 8 bins, x-coordinate")
plotbary(ph8,ptext=0.0023,coordi=2,col=col,nodesymbo=19,symbo="L",collines=col)
title(sub="a.2) 8 bins, y-coordinate")
col<-colobary(ph13$parent,paletti)
plotbary(ph13,ptext=0.005,coordi=1,col=col,nodesymbo=19,collines=col)
title(sub="b.1) 13 bins, x-coordinate")
plotbary(ph13,ptext=0.005,coordi=2,col=col,nodesymbo=19,collines=col)
title(sub="b.2) 13 bins, y-coordinate")
dev.off()
h<-1.4
N<-c(64,64)
Q<-80
pr12<-profkern(dendat,h,N,Q)
keva<-kereva(dendat,h,N)
dk12<-drawkern(keva$value,keva$index,N,dendat,h)
h<-0.8
N<-c(64,64)
Q<-80
pr07<-profkern(dendat,h,N,Q)
keva<-kereva(dendat,h,N)
dk07<-drawkern(keva$value,keva$index,N,dendat,h)
h<-0.6
N<-c(64,64)
Q<-80
pr05<-profkern(dendat,h,N,Q)
keva<-kereva(dendat,h,N)
dk05<-drawkern(keva$value,keva$index,N,dendat,h)
postscript(file="~/artlst/kerletree.ps",horizontal=FALSE)
par(mfcol=c(3,2)) #par(cex.sub=1.4)
parold<-par(no.readonly=FALSE)$plt
par(plt=c(0.08,0.95,0.11,0.91))
persp(dk12$x,dk12$y,dk12$z,theta=130,phi=20,col="white",
xlab="x",ylab="y",zlab="",ticktype="detailed")
text(-0.3,-0.45,"a)",cex=2)
#title(sub="a.1) h=1.4")
persp(dk07$x,dk07$y,dk07$z,theta=130,phi=20,col="white",#box=F,
xlab="x",ylab="y",zlab="",ticktype="detailed")
text(-0.3,-0.45,"b)",cex=2)
#title(sub="b.1) h=0.8")
persp(dk05$x,dk05$y,dk05$z,theta=130,phi=20,col="white",#box=F,
xlab="x",ylab="y",zlab="",ticktype="detailed")
text(-0.3,-0.45,"c)",cex=2)
#title(sub="c.1) h=0.6")
par(plt=parold) #par(plt=c( 0.2, 0.85, 0.2, 0.7))
plotvolu(pr12,ptext=0.0015,symbo="K")
#title(sub="a.2) h=1.4")
plotvolu(pr07,ptext=0.0025,symbo="L")
#title(sub="b.2) h=0.8")
plotvolu(pr05,ptext=0.0028)
#title(sub="c.2) h=0.6")
dev.off()
paletti<-c("red","blue","green",
"orange","navy","darkgreen",
"orchid","aquamarine","turquoise",
"pink","violet","magenta","chocolate","cyan",
colors()[50:100])
postscript(file="~/artlst/kerlebary.ps",horizontal=FALSE)
par(mfrow=c(3,2))
col<-colobary(pr12$parent,paletti)
ptext<-0.0016
plotbary(pr12,coordi=1,ptext=ptext,symbo="K",col=col,collines=col)
#,nodesymbo=19)
title(sub="a.1) h=1.4, x-coordinate")
plotbary(pr12,coordi=2,ptext=ptext,symbo="K",col=col,collines=col)
title(sub="a.2) h=1.4, y-coordinate")
col<-colobary(pr07$parent,paletti)
plotbary(pr07,coordi=1,ptext=0.0025,symbo="L",col=col,collines=col)
title(sub="b.1) h=0.8, x-coordinate")
plotbary(pr07,coordi=2,ptext=0.0025,symbo="L",col=col,collines=col)
title(sub="b.2) h=0.8, y-coordinate")
col<-colobary(pr05$parent,paletti)
plotbary(pr05,coordi=1,ptext=0.0025,symbo="M",col=col,collines=col)
title(sub="c.1) h=0.6, x-coordinate")
plotbary(pr05,coordi=2,ptext=0.0025,symbo="M",col=col,collines=col)
title(sub="c.2) h=0.6, y-coordinate")
dev.off()
- Data
- Histogram Estimate
- Projections
We have three dimensional data which is a random sample from
an equal mixture of 4 standard Gaussians.
The means of the members of the mixture lie in vertices of a tetrahedron.
dist<-3 # determine the distance between vertices of the tetrahedron
d<-3
moodi<-4
M<-matrix(0,moodi,d)
height<-sqrt(3)/2 # sqrt(3)/2 = 0.8660254
len<-1/(2*sqrt(3)) # 1/(2*sqrt(3)) = 0.2886751
kor<-sqrt(2/3) # sqrt(2/3) = 0.8164966
M[1,]<-dist*c(1/2,0,0) # ( 1.5, 0.0, 0.0)
M[2,]<-dist*c(-1/2,0,0) # (-1.5, 0.0, 0.0)
M[3,]<-dist*c(0,height,0) # ( 0.0, 2.6, 0.0)
M[4,]<-dist*c(0,len,kor) # ( 0.0, 0.9, 2.4)
sig<-matrix(1,moodi,d)
p0<-1/moodi
p<-p0*rep(1,moodi)
n<-5000
siemen<-16
dendat<-simmix(n,M,sig,p,siemen)
binlkm<-9
ph<-profhist(dendat,binlkm)
# save(file="/home/jsk/artlst/mix4c3dhis.var",list=c("ph"))
# load(file="/home/jsk/artlst/mix4c3dhis.var")
paletti<-c("red","blue","green",
"orange","navy","darkgreen",
"orchid","aquamarine","turquoise",
"pink","violet","magenta","chocolate","cyan",
colors()[50:100])
col<-colobary(ph$parent,paletti)
postscript(file="/home/jsk/artlst/mix4c3dhis.ps",horizontal=FALSE)
par(mfrow=c(3,2))
plotvolu(ph,modelabel=FALSE)
title(sub="a.1) volume plot")
#segments(170,0.009,200,0.009)
#segments(170,0.009,170,0.015)
#segments(200,0.009,200,0.015)
plotvolu(ph,ptext=0.00025,cutlev=0.009)
title(sub="a.2) upper levels of a.1)")
plotbary(ph,coordi=1,ptext=0.0006,col=col,collines=col)
title(sub="b.1) barycenter plot, 1st coordinate")
plotbary(ph,coordi=2,ptext=0.0006,col=col,collines=col)
title(sub="b.2) barycenter plot, 2nd coordinate")
plotbary(ph,coordi=3,ptext=0.0006,col=col,collines=col)
title(sub="b.3) barycenter plot, 3rd coordinate")
dev.off()
ph$center[,moodilkm(ph$parent)$modloc]
[,1] [,2] [,3] [,4]
[1,] -0.02209665 0.99589608 -1.04008939 -0.02209665
[2,] 1.07894346 0.04257649 0.04257649 3.15167739
[3,] 2.06719101 0.15173719 0.15173719 0.15173719
em<-excmas(ph)
plotprof(ph,modelabel=FALSE,length=em)
plotprof(ph,ptext=0.0002,cutlev=0.002,length=em)
Let us denote coordinates as (x,y,z).
dist<-3
dproj<-2
sig<-matrix(1,moodi,dproj)
Projection to x-y plane.
Mxy<-matrix(0,moodi,dproj)
Mxy[1,]<-dist*c(1/2,0)
Mxy[2,]<-dist*c(-1/2,0)
Mxy[3,]<-dist*c(0,height)
Mxy[4,]<-dist*c(0,len)
dmxy<-drawmix(Mxy,sig,p,plkm=60)
persp(dmxy$x,dmxy$y,dmxy$z,phi=30,theta=130,box=FALSE,col="white")
Projection to x-z plane.
Mxz<-matrix(0,moodi,dproj)
Mxz[1,]<-dist*c(1/2,0)
Mxz[2,]<-dist*c(-1/2,0)
Mxz[3,]<-dist*c(0,0)
Mxz[4,]<-dist*c(0,kor)
dmxz<-drawmix(Mxz,sig,p,plkm=60)
persp(dmxz$x,dmxz$y,dmxz$z,phi=30,theta=130,box=FALSE,col="white")
Projection to y-z plane.
Myz<-matrix(0,moodi,dproj)
Myz[1,]<-dist*c(0,0)
Myz[2,]<-dist*c(0,0)
Myz[3,]<-dist*c(height,0)
Myz[4,]<-dist*c(len,kor)
dmyz<-drawmix(Myz,sig,p,plkm=60)
persp(dmyz$x,dmyz$y,dmyz$z,phi=30,theta=130,box=FALSE,col="white")
postscript("/home/jsk/artlst/proj3d.ps",horizontal=FALSE)
par(mfrow=c(2,2))
par(plt=c(0,1,0,1))
persp(dmxy$x,dmxy$y,dmxy$z,phi=30,theta=130,box=FALSE,col="white")
#title(sub="a) projection of density to x-y plane, D=3")
text(-0.3,-0.4,"a)",cex=1.5)
persp(dmxz$x,dmxz$y,dmxz$z,phi=30,theta=130,box=FALSE,col="white")
#title(sub="b) projection of density to x-z plane, D=3")
text(-0.3,-0.4,"b)",cex=1.5)
persp(dmyz$x,dmyz$y,dmyz$z,phi=30,theta=130,box=FALSE,col="white")
#title(sub="c) projection of density to y-z plane, D=3")
text(-0.3,-0.4,"c)",cex=1.5)
dist<-4
Mxy<-matrix(0,moodi,dproj)
Mxy[1,]<-dist*c(1/2,0)
Mxy[2,]<-dist*c(-1/2,0)
Mxy[3,]<-dist*c(0,height)
Mxy[4,]<-dist*c(0,len)
dm<-drawmix(Mxy,sig,p,plkm=60)
persp(dm$x,dm$y,dm$z,phi=30,theta=130,box=FALSE,col="white")
#title(sub="d) projection of density to x-y plane, D=4")
text(-0.3,-0.4,"d)",cex=1.5)
dev.off()
- Data
- Kernel Estimate
- Projections
We have four dimensional data which is a random sample from
an equal mixture of 5 standard Gaussians.
The means of the members of the mixture lie in vertices of a pentahedron.
dist<-4 # determine the distance between vertices of the pentahedron
d<-4
moodi<-5
M<-matrix(0,moodi,d)
M[1,]<-dist*c(1/2, 0,0,0)
M[2,]<-dist*c(-1/2,0,0,0)
M[3,]<-dist*c(0,sqrt(3)/2,0,0)
M[4,]<-dist*c(0,1/(2*sqrt(3)),sqrt(2/3),0)
M[5,]<-dist*c(0,1/(2*sqrt(3)),1/(2*sqrt(6)),sqrt(15/24))
sig<-matrix(1,moodi,d)
p0<-1/moodi
p<-p0*rep(1,moodi)
n<-2000
siemen<-13
dendat<-simmix(n,M,sig,p,siemen)
h<-1
N<-c(16,16,16,16) #c(32,32,32,32)
Q<-40
pk<-profkern(dendat,h,N,Q) #,kernel="gauss")
# save(file="~/artlst/mix5c4dker.var",list=c("pk"))
# load(file="~/artlst/mix5c4dker.var")
paletti<-c("red","blue","green",
"orange","navy","darkgreen",
"orchid","aquamarine","turquoise",
"pink","violet","magenta","chocolate","cyan",
colors()[50:100])
col<-colobary(pk$parent,paletti)
postscript(file="/home/jsk/artlst/mix5c4dker.ps",horizontal=FALSE)
par(mfrow=c(3,2))
plotvolu(pk,modelabel=FALSE)
title(sub="a.1) volume plot")
plotvolu(pk,ptext=0.00021,cutlev=0.0017)
title(sub="a.2) upper levels of a.1)")
plotbary(pk,coordi=1,ptext=0.00017,col=col,collines=col)
title(sub="b.1) barycenter plot, 1st coordinate")
plotbary(pk,coordi=2,ptext=0.00017,col=col,collines=col)
title(sub="b.2) barycenter plot, 2nd coordinate")
plotbary(pk,coordi=3,ptext=0.00017,col=col,collines=col)
title(sub="b.3) barycenter plot, 3rd coordinate")
plotbary(pk,coordi=4,ptext=0.00017,col=col,collines=col)
title(sub="b.4) barycenter plot, 4th coordinate")
dev.off()
pk$center[,moodilkm(pk$parent)$modloc]
[,1] [,2] [,3] [,4] [,5]
[1,] 0.1786039 2.1518685 0.1786039 0.1786039 -1.7946607
[2,] 0.6909826 0.3563801 1.3601875 3.3678024 0.0217776
[3,] 0.9168697 0.1623964 3.1802897 -0.5920769 0.1623964
[4,] 3.0843982 -0.2972763 -0.2972763 -0.2972763 0.3790586
em<-excmas(pk)
plotprof(pk,modelabel=FALSE,length=em)
plotprof(pk,ptext=0.0002,cutlev=0.002,length=em)
Let us denote coordinates as (x,y,z,u).
dist<-4
dproj<-2
sig<-matrix(1,moodi,dproj)
Projection to x-y plane.
M<-matrix(0,moodi,dproj)
M[1,]<-dist*c(1/2,0)
M[2,]<-dist*c(-1/2,0)
M[3,]<-dist*c(0,sqrt(3)/2)
M[4,]<-dist*c(0,1/(2*sqrt(3)))
M[5,]<-dist*c(0,1/(2*sqrt(3)))
dmxy<-drawmix(M,sig,p,plkm=60)
persp(dmxy$x,dmxy$y,dmxy$z,phi=30,theta=130,box=FALSE,col="white")
Projection to x-z plane.
M<-matrix(0,moodi,dproj)
M[1,]<-dist*c(1/2,0)
M[2,]<-dist*c(-1/2,0)
M[3,]<-dist*c(0,0)
M[4,]<-dist*c(0,sqrt(2/3))
M[5,]<-dist*c(0,1/(2*sqrt(6)))
dmxz<-drawmix(M,sig,p,plkm=60)
persp(dmxz$x,dmxz$y,dmxz$z,phi=30,theta=130,box=FALSE,col="white")
Projection to x-u plane.
M<-matrix(0,moodi,dproj)
M[1,]<-dist*c(1/2,0)
M[2,]<-dist*c(-1/2,0)
M[3,]<-dist*c(0,0)
M[4,]<-dist*c(0,0)
M[5,]<-dist*c(0,sqrt(15/24))
dmxu<-drawmix(M,sig,p,plkm=60)
persp(dmxu$x,dmxu$y,dmxu$z,phi=30,theta=130,box=FALSE,col="white")
Projection to y-z plane.
M<-matrix(0,moodi,dproj)
M[1,]<-dist*c(0,0)
M[2,]<-dist*c(0,0)
M[3,]<-dist*c(sqrt(3)/2,0)
M[4,]<-dist*c(1/(2*sqrt(3)),sqrt(2/3))
M[5,]<-dist*c(1/(2*sqrt(3)),1/(2*sqrt(6)))
dmyz<-drawmix(M,sig,p,plkm=60)
persp(dmyz$x,dmyz$y,dmyz$z,phi=30,theta=130,box=FALSE,col="white")
Projection to y-u plane.
M<-matrix(0,moodi,dproj)
M[1,]<-dist*c(0,0)
M[2,]<-dist*c(0,0)
M[3,]<-dist*c(sqrt(3)/2,0)
M[4,]<-dist*c(1/(2*sqrt(3)),0)
M[5,]<-dist*c(1/(2*sqrt(3)),sqrt(15/24))
dmyu<-drawmix(M,sig,p,plkm=60)
persp(dmyu$x,dmyu$y,dmyu$z,phi=30,theta=130,box=FALSE,col="white")
Projection to z-u plane.
M<-matrix(0,moodi,dproj)
M[1,]<-dist*c(0,0)
M[2,]<-dist*c(0,0)
M[3,]<-dist*c(0,0)
M[4,]<-dist*c(sqrt(2/3),0)
M[5,]<-dist*c(1/(2*sqrt(6)),sqrt(15/24))
dmzu<-drawmix(M,sig,p,plkm=60)
persp(dmzu$x,dmzu$y,dmzu$z,phi=30,theta=130,box=FALSE,col="white")
postscript("/home/jsk/artlst/proj4d.ps")#,horizontal=FALSE)
par(mfrow=c(2,3))
par(plt=c(0,1,0,1))
persp(dmxy$x,dmxy$y,dmxy$z,phi=30,theta=130,box=FALSE,col="white")
text(-0.3,-0.4,"a)",cex=1.5)
persp(dmxz$x,dmxz$y,dmxz$z,phi=30,theta=130,box=FALSE,col="white")
text(-0.3,-0.4,"b)",cex=1.5)
persp(dmxu$x,dmxu$y,dmxu$z,phi=30,theta=130,box=FALSE,col="white")
text(-0.3,-0.4,"c)",cex=1.5)
persp(dmyz$x,dmyz$y,dmyz$z,phi=30,theta=130,box=FALSE,col="white")
text(-0.3,-0.4,"d)",cex=1.5)
persp(dmyu$x,dmyu$y,dmyu$z,phi=30,theta=130,box=FALSE,col="white")
text(-0.3,-0.4,"e)",cex=1.5)
persp(dmzu$x,dmzu$y,dmzu$z,phi=30,theta=130,box=FALSE,col="white")
text(-0.3,-0.4,"f)",cex=1.5)
dev.off()
postscript(file="~/artlst/histopersp.ps") #,horizontal=FALSE)
par(mfcol=c(2,2))
parold<-par(no.readonly=FALSE)$plt
par(plt=c(0.08,0.95,0.11,0.91))
#par(plt=c(0.05,0.95,0.07,0.93))
persp(dh8$x,dh8$y,dh8$z,col="white",theta=130,phi=20,
xlab="x",ylab="y",zlab="",ticktype="detailed") #,shade=0.5)
text(-0.3,-0.5,"a)",cex=1.5)
#title(sub="a) 8 bins")
persp(dh13$x,dh13$y,dh13$z,col="white",theta=130,phi=20,
xlab="x",ylab="y",zlab="",ticktype="detailed")
text(-0.3,-0.5,"b)",cex=1.5)
#title(sub="b) 13 bins")
par(plt=parold)
#par(plt=c( 0.1172537, 0.9399432, 0.1457143, 0.8828571))
plotvolu(ph8,ptext=0.0025,symbo="L")
plotvolu(ph13,ptext=0.004)
dev.off()
paletti<-c("red","blue","green",
"orange","navy","darkgreen",
"orchid","aquamarine","turquoise",
"pink","violet","magenta","chocolate","cyan",
colors()[50:100])
postscript(file="~/artlst/histobary.ps") #,horizontal=FALSE)
par(mfrow=c(2,2))
col<-colobary(ph8$parent,paletti)
plotbary(ph8,ptext=0.0023,coordi=1,col=col,nodesymbo=19,symbo="L",collines=col)
title(sub="a.1) 8 bins, x-coordinate")
plotbary(ph8,ptext=0.0023,coordi=2,col=col,nodesymbo=19,symbo="L",collines=col)
title(sub="a.2) 8 bins, y-coordinate")
col<-colobary(ph13$parent,paletti)
plotbary(ph13,ptext=0.005,coordi=1,col=col,nodesymbo=19,collines=col)
title(sub="b.1) 13 bins, x-coordinate")
plotbary(ph13,ptext=0.005,coordi=2,col=col,nodesymbo=19,collines=col)
title(sub="b.2) 13 bins, y-coordinate")
dev.off()
