• QR Code (size=2)
• QR Code (size=3)
• QR Code (size=4)
• QR Code (size=5)
• QR Code (size=6)

9356
9337
10149
9788
9770
9911
9429
8775
10189
10529
9914
9790
9625
9729
10589
9611
9388
9510
9690
8434
9844
10601
9942
10229
9381
9635
11228
9999
10089
11622
10533
9965
11567
11321
11686
11747
10595
10751
12199
11690
10978
11753
10839
10518
12183
11967
12363
12359
12162
12096
14325
12670
13865
13563
12734
12464
13389
13961
14088
13143
13413
13579
15388
13708
14689
14883
13991
13854
14364
15672
15904
14016

14271
14013
15912
14290
14744
14721
13918
13263
15660
15629
15113
14526
15132
14908
16167
14122
13985
14236
13921
12394
15454
16146
15107
14593
14695
14513
17071
15179
15460
17173
15938
15003
18216
17847
18029
17281
16706
16750
18912
17763
16736
18061
16713
16769
19514
19251
19951
19052
19555
19083
22534
18854
19801
20346
18169
19087
20842
21602
22360
20334
21215
20530
23152
20134
21193
21628
20823
20493
22106
24178
24958
21620

par1 <- as.numeric(par1) #factor
if (par2 == 'TRUE') par2 <- TRUE else par2 <- FALSE
if (par3 == 'TRUE') par3 <- TRUE else par3 <- FALSE
library(rpart)
compute.bagplot<-function(x,y,
factor=3, # expanding factor for bag to get the loop
approx.limit=300, # limit
dkmethod=2, # in 1:2; there are two methods for approximating the bag
precision=1, # controls precisionn of computation
verbose=FALSE,debug.plots='no' # tools for debugging
){
win<-function(dx,dy){  atan2(y=dy,x=dx) }
out.of.polygon<-function(xy,pg){
if(nrow(pg)==1) return(pg)
pgcenter<-apply(pg,2,mean)
pg<-cbind(pg[,1]-pgcenter[1],pg[,2]-pgcenter[2])
xy<-cbind(xy[,1]-pgcenter[1],xy[,2]-pgcenter[2])
extr<-rep(FALSE,nrow(xy))
for(i in seq(nrow(xy))){
alpha<-sort((win(xy[i,1]-pg[,1],xy[i,2]-pg[,2]))%%(2*pi))
extr[i]<-pipi<(alpha[1]+2*pi-alpha[length(alpha)])
}
extr
}
cut.z.pg<-function(zx,zy,p1x,p1y,p2x,p2y){
a2<-(p2y-p1y)/(p2x-p1x); a1<-zy/zx
sx<-(p1y-a2*p1x)/(a1-a2); sy<-a1*sx
sxy<-cbind(sx,sy)
h<-any(is.nan(sxy))||any(is.na(sxy))||any(Inf==abs(sxy))
if(h){
if(!exists('verbose')) verbose<-FALSE
if(verbose) cat('special')
h<-0==(a1-a2) & sign(zx)==sign(p1x)
sx<-ifelse(h,p1x,sx); sy<-ifelse(h,p1y,sy)
h<-0==(a1-a2) & sign(zx)!=sign(p1x)
sx<-ifelse(h,p2x,sx); sy<-ifelse(h,p2y,sy)
h<-p1x==p2x & zx!=p1x & p1x!=0
sx<-ifelse(h,p1x,sx); sy<-ifelse(h,zy*p1x/zx,sy)
h<-p1x==p2x & zx!=p1x & p1x==0
sx<-ifelse(h,p1x,sx); sy<-ifelse(h,0,sy)
h<-p1x==p2x & zx==p1x & p1x==0 & sign(zy)==sign(p1y)
sx<-ifelse(h,p1x,sx); sy<-ifelse(h,p1y,sy)
h<-p1x==p2x & zx==p1x & p1x==0 & sign(zy)!=sign(p1y)
sx<-ifelse(h,p1x,sx); sy<-ifelse(h,p2y,sy)
h<-zx==p1x & zy==p1y; sx<-ifelse(h,p1x,sx); sy<-ifelse(h,p1y,sy)
h<-zx==p2x & zy==p2y; sx<-ifelse(h,p2x,sx); sy<-ifelse(h,p2y,sy)
h<-zx==0 & zy==0; sx<-ifelse(h,0,sx); sy<-ifelse(h,0,sy)
sxy<-cbind(sx,sy)
} # end of special cases
if(!exists('debug.plots')) debug.plots<-'no'
if(debug.plots=='all'){
segments(sxy[,1],sxy[,2],zx,zy,col='red')
segments(0,0,sxy[,1],sxy[,2],type='l',col='green',lty=2)
points(sxy,col='red')
}
return(sxy)
}
find.cut.z.pg<-function(z,pg,center=c(0,0),debug.plots='no'){
if(!is.matrix(z)) z<-rbind(z)
if(1==nrow(pg)) return(matrix(center,nrow(z),2,TRUE))
n.pg<-nrow(pg); n.z<-nrow(z)
z<-cbind(z[,1]-center[1],z[,2]-center[2])
pgo<-pg; pg<-cbind(pg[,1]-center[1],pg[,2]-center[2])
if(!exists('debug.plots')) debug.plots<-'no'
if(debug.plots=='all'){plot(rbind(z,pg,0),bty='n'); points(z,pch='p')
lines(c(pg[,1],pg[1,1]),c(pg[,2],pg[1,2]))}
apg<-win(pg[,1],pg[,2])
apg[is.nan(apg)]<-0; a<-order(apg); apg<-apg[a]; pg<-pg[a,]
az<-win(z[,1],z[,2])
segm.no<-apply((outer(apg,az,'<')),2,sum)
segm.no<-ifelse(segm.no==0,n.pg,segm.no)
next.no<-1+(segm.no %% length(apg))
cuts<-cut.z.pg(z[,1],z[,2],pg[segm.no,1],pg[segm.no,2],
pg[next.no,1],pg[next.no,2])
cuts<-cbind(cuts[,1]+center[1],cuts[,2]+center[2])
return(cuts)
}
hdepth.of.points<-function(tp,n){
n.tp<-nrow(tp)
tphdepth<-rep(0,n.tp); dpi<-2*pi-0.000001
minusplus<-c(rep(-1,n),rep(1,n))
for(j in 1:n.tp) {
dx<-tp[j,1]-xy[,1]; dy<-tp[j,2]-xy[,2]
a<-win(dx,dy)+pi; h<-a<10;a<-a[h]; ident<-sum(!h)
init<-sum(a < pi); a.shift<-(a+pi) %% dpi
h<-cumsum(minusplus[order(c(a,a.shift))])
tphdepth[j]<-init+min(h)+1
}
tphdepth
}
expand.hull<-function(pg,k){
resolution<-floor(20*precision)
pg0<-xy[hdepth==1,]
pg0<-pg0[chull(pg0[,1],pg0[,2]),]
end.points<-find.cut.z.pg(pg,pg0,center=center,debug.plots=debug.plots)
lam<-((0:resolution)^1)/resolution^1
pg.new<-pg
for(i in 1:nrow(pg)){
tp<-cbind(pg[i,1]+lam*(end.points[i,1]-pg[i,1]),
pg[i,2]+lam*(end.points[i,2]-pg[i,2]))
hd.tp<-hdepth.of.points(tp,nrow(xy))
ind<-max(sum(hd.tp>=k),1)
if(indk &&
tp<-cbind(tp[ind,1]+lam*(tp[ind+1,1]-tp[ind,1]),
tp[ind,2]+lam*(tp[ind+1,2]-tp[ind,2]))
hd.tp<-hdepth.of.points(tp,nrow(xy))
ind<-max(sum(hd.tp>=k),1)
}
pg.new[i,]<-tp[ind,]
}
pg.new<-pg.new[chull(pg.new[,1],pg.new[,2]),]
pg.add<-0.5*(pg.new+rbind(pg.new[-1,],pg.new[1,]))
end.points<-find.cut.z.pg(pg,pg0,center=center)
for(i in 1:nrow(pg.add)){
tp<-cbind(pg.add[i,1]+lam*(end.points[i,1]-pg.add[i,1]),
pg.add[i,2]+lam*(end.points[i,2]-pg.add[i,2]))
hd.tp<-hdepth.of.points(tp,nrow(xy))
ind<-max(sum(hd.tp>=k),1)
if(indk &&
tp<-cbind(tp[ind,1]+lam*(tp[ind+1,1]-tp[ind,1]),
tp[ind,2]+lam*(tp[ind+1,2]-tp[ind,2]))
hd.tp<-hdepth.of.points(tp,nrow(xy))
ind<-max(sum(hd.tp>=k),1)
}
pg.add[i,]<-tp[ind,]
}
pg.new<-rbind(pg.new,pg.add)
pg.new<-pg.new[chull(pg.new[,1],pg.new[,2]),]
}
cut.p.sl.p.sl<-function(xy1,m1,xy2,m2){
sx<-(xy2[2]-m2*xy2[1]-xy1[2]+m1*xy1[1])/(m1-m2)
sy<-xy1[2]-m1*xy1[1]+m1*sx
if(!is.nan(sy)) return( c(sx,sy) )
if(abs(m1)==Inf) return( c(xy1[1],xy2[2]+m2*(xy1[1]-xy2[1])) )
if(abs(m2)==Inf) return( c(xy2[1],xy1[2]+m1*(xy2[1]-xy1[1])) )
}
pos.to.pg<-function(z,pg,reverse=FALSE){
if(reverse){
int.no<-apply(outer(pg[,1],z[,1],'>='),2,sum)
zy.on.pg<-pg[int.no,2]+pg[int.no,3]*(z[,1]-pg[int.no,1])
}else{
int.no<-apply(outer(pg[,1],z[,1],'<='),2,sum)
zy.on.pg<-pg[int.no,2]+pg[int.no,3]*(z[,1]-pg[int.no,1])
}
ifelse(z[,2]}
xydata<-if(missing(y)) x else cbind(x,y)
if(is.data.frame(xydata)) xydata<-as.matrix(xydata)
very.large.data.set<-nrow(xydata)>approx.limit
if(!exists('.Random.seed')) set.seed(13)
save.seed<-.Random.seed
if(very.large.data.set){
ind<-sample(seq(nrow(xydata)),size=approx.limit)
xy<-xydata[ind,]
} else xy<-xydata
n<-nrow(xy)
points.in.bag<-floor(n/2)
assign('.Random.seed',save.seed,env=.GlobalEnv)
if(verbose) cat('end of initialization')
prdata<-prcomp(xydata)
is.one.dim<-(min(prdata[[1]])/max(prdata[[1]]))<0.0001
if(is.one.dim){
if(verbose) cat('data set one dimensional')
center<-colMeans(xydata)
res<-list(xy=xy,xydata=xydata,prdata=prdata,is.one.dim=is.one.dim,center=center)
class(res)<-'bagplot'
return(res)
}
if(verbose) cat('data not linear')
dx<-(outer(xy[,1],xy[,1],'-'))
dy<-(outer(xy[,2],xy[,2],'-'))
alpha<-atan2(y=dy,x=dx); diag(alpha)<-1200
for(j in 1:n) alpha[,j]<-sort(alpha[,j])
alpha<-alpha[-n,] ; m<-n-1
if(debug.plots=='all'){
plot(xy,bty='n'); xdelta<-abs(diff(range(xy[,1]))); dx<-xdelta*.3
for(j in 1:n) {
p<-xy[j,]; dy<-dx*tan(alpha[,j])
segments(p[1]-dx,p[2]-dy,p[1]+dx,p[2]+dy,col=j)
text(p[1]-xdelta*.02,p[2],j,col=j)
}
}
if(verbose) print('end of computation of angles')
hdepth<-rep(0,n); dpi<-2*pi-0.000001
minusplus<-c(rep(-1,m),rep(1,m))
for(j in 1:n) {
a<-alpha[,j]+pi; h<-a<10; a<-a[h]; init<-sum(a < pi) # hallo
a.shift<-(a+pi) %% dpi
h<-cumsum(minusplus[order(c(a,a.shift))])
hdepth[j]<-init+min(h)+1 #  or do we have to count identical points?:
}
if(verbose){print('end of computation of hdepth:'); print(hdepth)}
if(debug.plots=='all'){
plot(xy,bty='n')
xdelta<-abs(diff(range(xy[,1]))); dx<-xdelta*.1
for(j in 1:n) {
a<-alpha[,j]+pi; a<-a[a<10]; init<-sum(a < pi)
a.shift<-(a+pi) %% dpi
h<-cumsum(minusplus[ao<-(order(c(a,a.shift)))])
no<-which((init+min(h)) == (init+h))[1]
p<-xy[j,]; dy<-dx*tan(alpha[,j])
segments(p[1]-dx,p[2]-dy,p[1]+dx,p[2]+dy,col=j,lty=3)
dy<-dx*tan(c(sort(a),sort(a))[no])
segments(p[1]-5*dx,p[2]-5*dy,p[1]+5*dx,p[2]+5*dy,col='black')
text(p[1]-xdelta*.02,p[2],hdepth[j],col=1,cex=2.5)
}
}
hd.table<-table(sort(hdepth))
d.k<-cbind(dk=rev(cumsum(rev(hd.table))),
k =as.numeric(names(hd.table)))
k.1<-sum(points.in.bagif(nrow(d.k)>1){
k<-d.k[k.1+1,2]
} else {
k<-d.k[k.1,2]
}
if(verbose){cat('counts of members of dk:'); print(hd.table)}
if(verbose){cat('end of computation of k, k=',k)}
center<-apply(xy[which(hdepth==max(hdepth)),,drop=FALSE],2,mean)
hull.center<-NULL
if(102){
n.p<-floor(c(32,16,8)[1+(n>50)+(n>200)]*precision)
cands<-xy[rev(order(hdepth))[1:6],]
cands<-cands[chull(cands[,1],cands[,2]),]; n.c<-nrow(cands)
xyextr<-rbind(apply(cands,2,min),apply(cands,2,max))
h1<-seq(xyextr[1,1],xyextr[2,1],length=n.p)
h2<-seq(xyextr[1,2],xyextr[2,2],length=n.p)
tp<-cbind(matrix(h1,n.p,n.p)[1:n.p^2],
matrix(h2,n.p,n.p,TRUE)[1:n.p^2])
tphdepth<-hdepth.of.points(tp,n)
hull.center<-tp[which(tphdepth>=(max(tphdepth))),,drop=FALSE]
center<-apply(hull.center,2,mean)
cands<-hull.center[chull(hull.center[,1],hull.center[,2]),,drop=F]
xyextr<-rbind(apply(cands,2,min),apply(cands,2,max))
xydel<-(xyextr[2,]-xyextr[1,])/n.p
xyextr<-rbind(xyextr[1,]-xydel,xyextr[2,]+xydel)
h1<-seq(xyextr[1,1],xyextr[2,1],length=n.p)
h2<-seq(xyextr[1,2],xyextr[2,2],length=n.p)
tp<-cbind(matrix(h1,n.p,n.p)[1:n.p^2],
matrix(h2,n.p,n.p,TRUE)[1:n.p^2])
tphdepth<-hdepth.of.points(tp,n)
hull.center<-tp[which(tphdepth>=max(tphdepth)),,drop=FALSE]
center<-apply(hull.center,2,mean)
hull.center<-hull.center[chull(hull.center[,1],hull.center[,2]),]
if(verbose){cat('hull.center',hull.center); print(table(tphdepth)) }
}
if(verbose) cat('center depth:',hdepth.of.points(rbind(center),n))
if(verbose){print('end of computation of center'); print(center)}
if(dkmethod==1){
xyi<-xy[hdepth>=k,,drop=FALSE]
pdk<-xyi[chull(xyi[,1],xyi[,2]),,drop=FALSE]
xyo<-xy[hdepth>=(k-1),,drop=FALSE]
pdk.1<-xyo[chull(xyo[,1],xyo[,2]),,drop=FALSE]
if(verbose)cat('hull computed:')
if(debug.plots=='all'){
plot(xy,bty='n')
h<-rbind(pdk,pdk[1,]); lines(h,col='red',lty=2)
h<-rbind(pdk.1,pdk.1[1,]);lines(h,col='blue',lty=3)
points(center[1],center[2],pch=8,col='red')
}
exp.dk<-expand.hull(pdk,k)
exp.dk.1<-expand.hull(exp.dk,k-1) # pdk.1,k-1,20)
}else{
num<-floor(c(417,351,171,85,67,43)[sum(n>c(1,50,100,150,200,250))]*precision)
num.h<-floor(num/2); angles<-seq(0,pi,length=num.h)
ang<-tan(pi/2-angles)
xym<-apply(xy,2,mean); xysd<-apply(xy,2,sd)
xyxy<-cbind((xy[,1]-xym[1])/xysd[1],(xy[,2]-xym[2])/xysd[2])
kkk<-k
ia<-1; a<-angles[ia]; xyt<-xyxy%*%c(cos(a),-sin(a)); xyto<-order(xyt)
ind.k <-xyto[kkk]; cutp<-c(xyxy[ind.k,1],-10)
dxy<-diff(range(xyxy))
pg<-rbind(c(cutp[1],-dxy,Inf),c(cutp[1],dxy,NA))
ind.kk<-xyto[n+1-kkk]; cutpl<-c(xyxy[ind.kk,1],10)
pgl<-rbind(c(cutpl[1],dxy,Inf),c(cutpl[1],-dxy,NA))
if(debug.plots=='all'){ plot(xyxy,type='p',bty='n')
}
for(ia in seq(angles)[-1]){
a<-angles[ia]; angtan<-ang[ia]; xyt<-xyxy%*%c(cos(a),-sin(a)); xyto<-order(xyt)
ind.k <-xyto[kkk]; ind.kk<-xyto[n+1-kkk]; pnew<-xyxy[ind.k,]; pnewl<-xyxy[ind.kk,]
if(debug.plots=='all') points(pnew[1],pnew[2],col='red')
if(abs(angtan)>1e10){  ###  cat('y=c case')
pg.no<-sum(pg[,1]cutp<-c(pnew[1],pg[pg.no,2]+pg[pg.no,3]*(pnew[1]-pg[pg.no,1]))
pg.nol<-sum(pgl[,1]>=pnewl[1])
cutpl<-c(pnewl[1],pgl[pg.nol,2]+pgl[pg.nol,3]*(pnewl[1]-pgl[pg.nol,1]))
}else{    ### cat('normal case')
pg.inter<-pg[,2]-angtan*pg[,1]; pnew.inter<-pnew[2]-angtan*pnew[1]
pg.no<-sum(pg.intercutp<-cut.p.sl.p.sl(pnew,ang[ia],pg[pg.no,1:2],pg[pg.no,3])
pg.interl<-pgl[,2]-angtan*pgl[,1]; pnew.interl<-pnewl[2]-angtan*pnewl[1]
pg.nol<-sum(pg.interl>pnew.interl)
cutpl<-cut.p.sl.p.sl(pnewl,angtan,pgl[pg.nol,1:2],pgl[pg.nol,3])
}
pg<-rbind(pg[1:pg.no,],c(cutp,angtan),c(cutp[1]+dxy, cutp[2]+angtan*dxy,NA))
pgl<-rbind(pgl[1:pg.nol,],c(cutpl,angtan),c(cutpl[1]-dxy, cutpl[2]-angtan*dxy,NA))
if(debug.plots=='all'){
points(pnew[1],pnew[2],col='red')
hx<-xyxy[ind.k,c(1,1)]; hy<-xyxy[ind.k,c(2,2)]
segments(hx,hy,c(10,-10),hy+ang[ia]*(c(10,-10)-hx),lty=2)
points(cutpl[1],cutpl[2],col='red')
hx<-xyxy[ind.kk,c(1,1)]; hy<-xyxy[ind.kk,c(2,2)]
segments(hx,hy,c(10,-10),hy+ang[ia]*(c(10,-10)-hx),lty=2)
}
}
pg<-pg[-nrow(pg),][-1,,drop=F]; pgl<-pgl[-nrow(pgl),][-1,,drop=FALSE]
indl<-pos.to.pg(pgl,pg); indu<-pos.to.pg(pg,pgl,TRUE)
npg<-nrow(pg); npgl<-nrow(pgl)
rnuml<-rnumu<-lnuml<-lnumu<-0; sl<-pg[1,1:2]; sr<-pgl[1,1:2]
if(indl[1]=='higher'&indu[npg]=='lower'){
rnuml<-which(indl=='lower')[1]-1; xyl<-pgl[rnuml,] #
rnumu<-which(rev(indu=='higher'))[1]; xyu<-pg[npg+1-rnumu,] #
sr<-cut.p.sl.p.sl(xyl[1:2],xyl[3],xyu[1:2],xyu[3])
}
if(indl[npgl]=='higher'&indu[1]=='lower'){
lnuml<-which(rev(indl=='lower'))[1]; xyl<-pgl[npgl+1-lnuml,] #
lnumu<-which(indu=='higher')[1]-1; xyu<-pg[lnumu,] #?
sl<-cut.p.sl.p.sl(xyl[1:2],xyl[3],xyu[1:2],xyu[3])
}
pgl<-pgl[(rnuml+1):(npgl-lnuml),1:2,drop=FALSE]
pg <-pg [(lnumu+1):(npg -rnumu),1:2,drop=FALSE]
pg<-rbind(pg,sr,pgl,sl)
pg<-pg[chull(pg[,1],pg[,2]),]
if(debug.plots=='all') lines(rbind(pg,pg[1,]),col='red')
exp.dk<-cbind(pg[,1]*xysd[1]+xym[1],pg[,2]*xysd[2]+xym[2])
if(kkk>1) kkk<-kkk-1
ia<-1; a<-angles[ia]; xyt<-xyxy%*%c(cos(a),-sin(a)); xyto<-order(xyt)
ind.k <-xyto[kkk]; cutp<-c(xyxy[ind.k,1],-10)
dxy<-diff(range(xyxy))
pg<-rbind(c(cutp[1],-dxy,Inf),c(cutp[1],dxy,NA))
ind.kk<-xyto[n+1-kkk]; cutpl<-c(xyxy[ind.kk,1],10)
pgl<-rbind(c(cutpl[1],dxy,Inf),c(cutpl[1],-dxy,NA))
if(debug.plots=='all'){ plot(xyxy,type='p',bty='n')
}
for(ia in seq(angles)[-1]){
a<-angles[ia]; angtan<-ang[ia]; xyt<-xyxy%*%c(cos(a),-sin(a)); xyto<-order(xyt)
ind.k <-xyto[kkk]; ind.kk<-xyto[n+1-kkk]; pnew<-xyxy[ind.k,]; pnewl<-xyxy[ind.kk,]
if(debug.plots=='all') points(pnew[1],pnew[2],col='red')
if(abs(angtan)>1e10){  ###  cat('y=c case')
pg.no<-sum(pg[,1]cutp<-c(pnew[1],pg[pg.no,2]+pg[pg.no,3]*(pnew[1]-pg[pg.no,1]))
pg.nol<-sum(pgl[,1]>=pnewl[1])
cutpl<-c(pnewl[1],pgl[pg.nol,2]+pgl[pg.nol,3]*(pnewl[1]-pgl[pg.nol,1]))
}else{    ### cat('normal case')
pg.inter<-pg[,2]-angtan*pg[,1]; pnew.inter<-pnew[2]-angtan*pnew[1]
pg.no<-sum(pg.intercutp<-cut.p.sl.p.sl(pnew,ang[ia],pg[pg.no,1:2],pg[pg.no,3])
pg.interl<-pgl[,2]-angtan*pgl[,1]; pnew.interl<-pnewl[2]-angtan*pnewl[1]
pg.nol<-sum(pg.interl>pnew.interl)
cutpl<-cut.p.sl.p.sl(pnewl,angtan,pgl[pg.nol,1:2],pgl[pg.nol,3])
}
pg<-rbind(pg[1:pg.no,],c(cutp,angtan),c(cutp[1]+dxy, cutp[2]+angtan*dxy,NA))
pgl<-rbind(pgl[1:pg.nol,],c(cutpl,angtan),c(cutpl[1]-dxy, cutpl[2]-angtan*dxy,NA))
if(debug.plots=='all'){
points(pnew[1],pnew[2],col='red')
hx<-xyxy[ind.k,c(1,1)]; hy<-xyxy[ind.k,c(2,2)]
segments(hx,hy,c(10,-10),hy+ang[ia]*(c(10,-10)-hx),lty=2)
points(cutpl[1],cutpl[2],col='red')
hx<-xyxy[ind.kk,c(1,1)]; hy<-xyxy[ind.kk,c(2,2)]
segments(hx,hy,c(10,-10),hy+ang[ia]*(c(10,-10)-hx),lty=2)
}
}
pg<-pg[-nrow(pg),][-1,,drop=F]; pgl<-pgl[-nrow(pgl),][-1,,drop=FALSE]
indl<-pos.to.pg(pgl,pg); indu<-pos.to.pg(pg,pgl,TRUE)
npg<-nrow(pg); npgl<-nrow(pgl)
rnuml<-rnumu<-lnuml<-lnumu<-0; sl<-pg[1,1:2]; sr<-pgl[1,1:2]
if(indl[1]=='higher'&indu[npg]=='lower'){
rnuml<-which(indl=='lower')[1]-1; xyl<-pgl[rnuml,] #
rnumu<-which(rev(indu=='higher'))[1]; xyu<-pg[npg+1-rnumu,] #
sr<-cut.p.sl.p.sl(xyl[1:2],xyl[3],xyu[1:2],xyu[3])
}
if(indl[npgl]=='higher'&indu[1]=='lower'){
lnuml<-which(rev(indl=='lower'))[1]; xyl<-pgl[npgl+1-lnuml,] #
lnumu<-which(indu=='higher')[1]-1; xyu<-pg[lnumu,] #?
sl<-cut.p.sl.p.sl(xyl[1:2],xyl[3],xyu[1:2],xyu[3])
}
pgl<-pgl[(rnuml+1):(npgl-lnuml),1:2,drop=FALSE]
pg <-pg [(lnumu+1):(npg -rnumu),1:2,drop=FALSE]
pg<-rbind(pg,sr,pgl,sl)
pg<-pg[chull(pg[,1],pg[,2]),]
if(debug.plots=='all') lines(rbind(pg,pg[1,]),col='red')
exp.dk.1<-cbind(pg[,1]*xysd[1]+xym[1],pg[,2]*xysd[2]+xym[2])
}
lambda<-if(nrow(d.k)==1) 0.5 else
(n/2-d.k[k.1+1,1])/(d.k[k.1,1]-d.k[k.1+1,1])
if(verbose) cat('lambda',lambda)
cut.on.pdk.1<-find.cut.z.pg(exp.dk,  exp.dk.1,center=center)
cut.on.pdk  <-find.cut.z.pg(exp.dk.1,exp.dk,  center=center)
h1<-(1-lambda)*exp.dk+lambda*cut.on.pdk.1
h2<-(1-lambda)*cut.on.pdk+lambda*exp.dk.1
h<-rbind(h1,h2); hull.bag<-h[chull(h[,1],h[,2]),]
if(verbose)cat('bag completed:') #if(verbose)print(hull.bag)
if(debug.plots=='all'){   lines(hull.bag,col='red') }
hull.loop<-cbind(hull.bag[,1]-center[1],hull.bag[,2]-center[2])
hull.loop<-factor*hull.loop
hull.loop<-cbind(hull.loop[,1]+center[1],hull.loop[,2]+center[2])
if(verbose) cat('loop computed')
if(!very.large.data.set){
pxy.bag    <-xydata[hdepth>= k   ,,drop=FALSE]
pkt.cand   <-xydata[hdepth==(k-1),,drop=FALSE]
pkt.not.bag<-xydata[hdepth< (k-1),,drop=FALSE]
if(length(pkt.cand)>0){
outside<-out.of.polygon(pkt.cand,hull.bag)
if(sum(!outside)>0)
pxy.bag    <-rbind(pxy.bag,     pkt.cand[!outside,])
if(sum( outside)>0)
pkt.not.bag<-rbind(pkt.not.bag, pkt.cand[ outside,])
}
}else {
extr<-out.of.polygon(xydata,hull.bag)
pxy.bag    <-xydata[!extr,]
pkt.not.bag<-xydata[extr,,drop=FALSE]
}
if(length(pkt.not.bag)>0){
extr<-out.of.polygon(pkt.not.bag,hull.loop)
pxy.outlier<-pkt.not.bag[extr,,drop=FALSE]
if(0==length(pxy.outlier)) pxy.outlier<-NULL
pxy.outer<-pkt.not.bag[!extr,,drop=FALSE]
}else{
pxy.outer<-pxy.outlier<-NULL
}
if(verbose) cat('points of bag, outer points and outlier identified')
hull.loop<-rbind(pxy.outer,hull.bag)
hull.loop<-hull.loop[chull(hull.loop[,1],hull.loop[,2]),]
if(verbose) cat('end of computation of loop')
assign('.Random.seed',save.seed,env=.GlobalEnv)
res<-list(
center=center,
pxy.bag=pxy.bag,
pxy.outer=if(length(pxy.outer)>0) pxy.outer else NULL,
pxy.outlier=if(length(pxy.outlier)>0) pxy.outlier else NULL,
hull.center=hull.center,
hull.bag=hull.bag,
hull.loop=hull.loop,
hdepths=hdepth,
is.one.dim=is.one.dim,
prdata=prdata,
xy=xy,xydata=xydata
)
if(verbose) res<-c(res,list(exp.dk=exp.dk,exp.dk.1=exp.dk.1,hdepth=hdepth))
class(res)<-'bagplot'
return(res)
}
plot.bagplot<-function(bagplot.obj,
show.outlier=TRUE,# if TRUE outlier are shown
show.whiskers=TRUE, # if TRUE whiskers are shown
show.looppoints=TRUE, # if TRUE points in loop are shown
show.bagpoints=TRUE, # if TRUE points in bag are shown
show.loophull=TRUE, # if TRUE loop is shown
show.baghull=TRUE, # if TRUE bag is shown
add=FALSE, # if TRUE graphical elements are added to actual plot
pch=16,cex=.4,..., # to define further parameters of plot
verbose=FALSE # tools for debugging
){
win<-function(dx,dy){  atan2(y=dy,x=dx) }
cut.z.pg<-function(zx,zy,p1x,p1y,p2x,p2y){
a2<-(p2y-p1y)/(p2x-p1x); a1<-zy/zx
sx<-(p1y-a2*p1x)/(a1-a2); sy<-a1*sx
sxy<-cbind(sx,sy)
h<-any(is.nan(sxy))||any(is.na(sxy))||any(Inf==abs(sxy))
if(h){
if(!exists('verbose')) verbose<-FALSE
if(verbose) cat('special')
h<-0==(a1-a2) & sign(zx)==sign(p1x)
sx<-ifelse(h,p1x,sx); sy<-ifelse(h,p1y,sy)
h<-0==(a1-a2) & sign(zx)!=sign(p1x)
sx<-ifelse(h,p2x,sx); sy<-ifelse(h,p2y,sy)
h<-p1x==p2x & zx!=p1x & p1x!=0
sx<-ifelse(h,p1x,sx); sy<-ifelse(h,zy*p1x/zx,sy)
h<-p1x==p2x & zx!=p1x & p1x==0
sx<-ifelse(h,p1x,sx); sy<-ifelse(h,0,sy)
h<-p1x==p2x & zx==p1x & p1x==0 & sign(zy)==sign(p1y)
sx<-ifelse(h,p1x,sx); sy<-ifelse(h,p1y,sy)
h<-p1x==p2x & zx==p1x & p1x==0 & sign(zy)!=sign(p1y)
sx<-ifelse(h,p1x,sx); sy<-ifelse(h,p2y,sy)
h<-zx==p1x & zy==p1y; sx<-ifelse(h,p1x,sx); sy<-ifelse(h,p1y,sy)
h<-zx==p2x & zy==p2y; sx<-ifelse(h,p2x,sx); sy<-ifelse(h,p2y,sy)
h<-zx==0 & zy==0; sx<-ifelse(h,0,sx); sy<-ifelse(h,0,sy)
sxy<-cbind(sx,sy)
} # end of special cases
if(!exists('debug.plots')) debug.plots<-'no'
if(debug.plots=='all'){
segments(sxy[,1],sxy[,2],zx,zy,col='red')
segments(0,0,sxy[,1],sxy[,2],type='l',col='green',lty=2)
points(sxy,col='red')
}
return(sxy)
}
find.cut.z.pg<-function(z,pg,center=c(0,0),debug.plots='no'){
if(!is.matrix(z)) z<-rbind(z)
if(1==nrow(pg)) return(matrix(center,nrow(z),2,TRUE))
n.pg<-nrow(pg); n.z<-nrow(z)
z<-cbind(z[,1]-center[1],z[,2]-center[2])
pgo<-pg; pg<-cbind(pg[,1]-center[1],pg[,2]-center[2])
if(!exists('debug.plots')) debug.plots<-'no'
if(debug.plots=='all'){plot(rbind(z,pg,0),bty='n'); points(z,pch='p')
lines(c(pg[,1],pg[1,1]),c(pg[,2],pg[1,2]))}
apg<-win(pg[,1],pg[,2])
apg[is.nan(apg)]<-0; a<-order(apg); apg<-apg[a]; pg<-pg[a,]
az<-win(z[,1],z[,2])
segm.no<-apply((outer(apg,az,'<')),2,sum)
segm.no<-ifelse(segm.no==0,n.pg,segm.no)
next.no<-1+(segm.no %% length(apg))
cuts<-cut.z.pg(z[,1],z[,2],pg[segm.no,1],pg[segm.no,2],
pg[next.no,1],pg[next.no,2])
cuts<-cbind(cuts[,1]+center[1],cuts[,2]+center[2])
return(cuts)
}
for(i in seq(along=bagplot.obj))
eval(parse(text=paste(names(bagplot.obj)[i],'<-bagplot.obj[[',i,']]')))
if(is.one.dim){
if(verbose) cat('data set one dimensional')
prdata<-prdata[[2]];
trdata<-xydata%*%prdata; ytr<-mean(trdata[,2])
boxplotres<-boxplot(trdata[,1],plot=FALSE)
dy<-0.1*diff(range(stats<-boxplotres$stats))
dy<-0.05*mean(c(diff(range(xydata[,1])),
diff(range(xydata[,2]))))
segtr<-rbind(cbind(stats[2:4],ytr-dy,stats[2:4],ytr+dy),
cbind(stats[c(2,2)],ytr+c(dy,-dy),
stats[c(4,4)],ytr+c(dy,-dy)),
cbind(stats[c(2,4)],ytr,stats[c(1,5)],ytr))
segm<-cbind(segtr[,1:2]%*%t(prdata),
segtr[,3:4]%*%t(prdata))
if(!add) plot(xydata,type='n',bty='n',pch=16,cex=.2,...)
extr<-c(min(segm[6,3],segm[7,3]),max(segm[6,3],segm[7,3]))
extr<-extr+c(-1,1)*0.000001*diff(extr)
xydata<-xydata[xydata[,1]xydata[,1]>extr[2],,drop=FALSE]
if(0segments(segm[,1],segm[,2],segm[,3],segm[,4],)
return('one dimensional boxplot plottet')
}
if(!add) plot(xydata,type='n',pch=pch,cex=cex,bty='n',...)
if(verbose) text(xy[,1],xy[,2],paste(as.character(hdepth)),cex=2)
if(show.loophull){ # fill loop
h<-rbind(hull.loop,hull.loop[1,]); lines(h[,1],h[,2],lty=1)
polygon(hull.loop[,1],hull.loop[,2],col='#aaccff')
}
if(show.looppoints && length(pxy.outer)>0){ # points in loop
points(pxy.outer[,1],pxy.outer[,2],col='#3355ff',pch=pch,cex=cex)
}
if(show.baghull){ # fill bag
h<-rbind(hull.bag,hull.bag[1,]); lines(h[,1],h[,2],lty=1)
polygon(hull.bag[,1],hull.bag[,2],col='#7799ff')
}
if(show.bagpoints && length(pxy.bag)>0){ # points in bag
points(pxy.bag[,1],pxy.bag[,2],col='#000088',pch=pch,cex=cex)
}
if(show.whiskers && length(pxy.outer)>0){
debug.plots<-'not'
pkt.cut<-find.cut.z.pg(pxy.outer,hull.bag,center=center)
segments(pxy.outer[,1],pxy.outer[,2],pkt.cut[,1],pkt.cut[,2],col='red')
}
if(show.outlier && length(pxy.outlier)>0){ # points in loop
points(pxy.outlier[,1],pxy.outlier[,2],col='red',pch=pch,cex=cex)
}
if(exists('hull.center')&&length(hull.center)>2){
h<-rbind(hull.center,hull.center[1,]); lines(h[,1],h[,2],lty=1)
polygon(hull.center[,1],hull.center[,2],col='orange')
}
points(center[1],center[2],pch=8,col='red')
if(verbose){
h<-rbind(exp.dk,exp.dk[1,]); lines(h,col='blue',lty=2)
h<-rbind(exp.dk.1,exp.dk.1[1,]); lines(h,col='black',lty=2)
if(exists('tphdepth'))
text(tp[,1],tp[,2],as.character(tphdepth),col='green')
text(xy[,1],xy[,2],paste(as.character(hdepth)),cex=2)
points(center[1],center[2],pch=8,col='red')
}
'bagplot plottet'
}
bagplot<-function(x,y,
factor=3, # expanding factor for bag to get the loop
approx.limit=300, # limit
show.outlier=TRUE,# if TRUE outlier are shown
show.whiskers=TRUE, # if TRUE whiskers are shown
show.looppoints=TRUE, # if TRUE points in loop are shown
show.bagpoints=TRUE, # if TRUE points in bag are shown
show.loophull=TRUE, # if TRUE loop is shown
show.baghull=TRUE, # if TRUE bag is shown
create.plot=TRUE, # if TRUE a plot is created
add=FALSE, # if TRUE graphical elements are added to actual plot
pch=16,cex=.4,..., # to define further parameters of plot
dkmethod=2, # in 1:2; there are two methods for approximating the bag
precision=1, # controls precisionn of computation
verbose=FALSE,debug.plots='' # tools for debugging
){
bo<-compute.bagplot(x=x,y=y,factor=factor,approx.limit=approx.limit,
dkmethod=dkmethod,precision=precision,
verbose=verbose,debug.plots=debug.plots)
if(create.plot){
plot(bo,
show.outlier=show.outlier,
show.whiskers=show.whiskers,
show.looppoints=show.looppoints,
show.bagpoints=show.bagpoints,
show.loophull=show.loophull,
show.baghull=show.baghull,
add=add,pch=pch,cex=cex,...,
verbose=verbose
)
}
}
bitmap(file='test1.png')
bagplot(x=x, y=y, verbose=F, factor=par1, show.outlier=par2, show.whiskers=par3, show.baghull=T, dkmethod=2, show.loophull=T, precision=1, xlab=xlab, ylab=ylab, main=main)
box()
dev.off()