Send output to:
Browser Blue - Charts White
Browser Black/White
CSV
Data X:
611 594 595 591 589 584 573 567 569 621 629 628 612 595 597 593 590 580 574 573 573 620 626 620 588 566 557 561 549 532 526 511 499 555 565 542 527 510 514 517 508 493 490 469 478 528 534 518 506 502 516 528 533 536 537 524 536 587 597 581
Data Y:
120.9 119.6 125.9 116.1 107.5 116.7 112.5 113 126.4 114.1 112.5 112.4 113.1 116.3 111.7 118.8 116.5 125.1 113.1 119.6 114.4 114 117.8 117 120.9 115 117.3 119.4 114.9 125.8 117.6 117.6 114.9 121.9 117 106.4 110.5 113.6 114.2 125.4 124.6 120.2 120.8 111.4 124.1 120.2 125.5 116 117 105.7 102 106.4 96.9 107.6 98.8 101.1 105.7 104.6 103.2 101.6
Sample Range:
(leave blank to include all observations)
From:
To:
Factor (defining the loop)
Show outliers
TRUE
TRUE
FALSE
Show whiskers
TRUE
TRUE
FALSE
Chart options
Title:
Label y-axis:
Label x-axis:
R Code
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]<-pi<max(diff(alpha)) | pi<(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(ind<length(hd.tp)){ # hd.tp[ind]>k && 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(ind<length(hd.tp)){ # hd.tp[ind]>k && 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]<zy.on.pg, 'lower','higher') } 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.bag<d.k[,1]) if(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(10<nrow(xy)&&length(hd.table)>2){ 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]<pnew[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.inter<pnew.inter) cutp<-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]<pnew[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.inter<pnew.inter) cutp<-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]<extr[1] | xydata[,1]>extr[2],,drop=FALSE] if(0<nrow(xydata))points(xydata[,1],xydata[,2],pch=pch,cex=cex) segments(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()
Compute
Summary of computational transaction
Raw Input
view raw input (R code)
Raw Output
view raw output of R engine
Computing time
0 seconds
R Server
Big Analytics Cloud Computing Center
Click here to blog (archive) this computation