Send output to:
Browser Blue - Charts White
Browser Black/White
CSV
Data X:
100 114.242116 115.8697864 108.0366226 100.3051882 98.1688708 97.5584944 105.4933876 104.8830112 105.4933876 114.0386572 88.40284842 96.64292981 113.733469 110.681587 111.190234 103.1536114 102.2380468 102.4415056 115.4628688 107.9348932 113.5300102 116.8870804 90.13224822 101.2207528 117.0905392 120.0406918 113.3265514 109.155646 107.121058 107.121058 119.9389624 112.1057986 114.3438454 119.5320448 94.60834181 105.2899288 118.311292 122.0752798 116.276704 106.5106816 111.698881 114.5473042 116.3784334 117.7009156 116.6836216 120.4476094 96.54120041 105.595117 117.0905392 115.6663276 105.7985758 95.93082401 94.09969481 94.81180061 106.5106816 95.62563581 99.7965412 104.4760936 83.82502543
Data Y:
100 100 93.5483871 88.17204301 89.24731183 91.39784946 92.47311828 91.39784946 88.17204301 87.09677419 84.94623656 92.47311828 93.5483871 93.5483871 91.39784946 90.32258065 91.39784946 93.5483871 93.5483871 92.47311828 91.39784946 89.24731183 86.02150538 88.17204301 87.09677419 87.09677419 86.02150538 84.94623656 84.94623656 86.02150538 86.02150538 84.94623656 86.02150538 82.79569892 77.41935484 80.64516129 78.49462366 75.2688172 75.2688172 75.2688172 77.41935484 78.49462366 76.34408602 73.11827957 68.8172043 65.59139785 69.89247312 82.79569892 84.94623656 80.64516129 74.19354839 70.96774194 74.19354839 82.79569892 86.02150538 86.02150538 82.79569892 78.49462366 79.56989247 87.09677419
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
1 seconds
R Server
Big Analytics Cloud Computing Center
Click here to blog (archive) this computation