Send output to:
Browser Blue - Charts White
Browser Black/White
CSV
Data X:
8.6 8.5 8.3 7.8 7.8 8 8.6 8.9 8.9 8.6 8.3 8.3 8.3 8.4 8.5 8.4 8.6 8.5 8.5 8.5 8.5 8.5 8.5 8.5 8.5 8.5 8.5 8.5 8.6 8.4 8.1 8 8 8 8 7.9 7.8 7.8 7.9 8.1 8 7.6 7.3 7 6.8 7 7.1 7.2 7.1 6.9 6.7 6.7 6.6 6.9 7.3 7.5 7.3 7.1 6.9 7.1
Data Y:
562325 560854 555332 543599 536662 542722 593530 610763 612613 611324 594167 595454 590865 589379 584428 573100 567456 569028 620735 628884 628232 612117 595404 597141 593408 590072 579799 574205 572775 572942 619567 625809 619916 587625 565742 557274 560576 548854 531673 525919 511038 498662 555362 564591 541657 527070 509846 514258 516922 507561 492622 490243 469357 477580 528379 533590 517945 506174 501866 516141
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