program asmrshis cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Version from February 3, 1996 c programmed by Joachim Engel c c program for average shifted multiresolution histogram c c Call: asmrshis dataset m xmin xmax lambda c e.g. asmrshis faithful.dat 100 1.2 5. 5 c c dataset set of univariate data (may be unordered) c m number of shifts c xmin c xmax algorithm applies to interval [xmin,xmax] c lambda finest resolution level c c default for boundaries: smallest and largest observation c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc implicit undefined (a-z) c integer m,nmax,n,nnb,nplot,ny,i,j,lambda character*30 dataset parameter(nplot=1000,nmax=1000) integer w(nmax) double precision t(nplot),fhad(nplot),fhad2(nplot) double precision dy(nmax),y(nmax),bb(nmax),s(nmax) double precision xmin,xmax,dt,bw,to,delta,xmax1 data ny/1/,to/0/ c call parm(dataset,m,xmin,xmax,lambda) c open(unit=50,file='mrshis.est') open(unit=60,file=dataset) do 100 i=1,nmax read(60,*,end=3) y(i) 100 continue 3 n=i-1 call bsort(y,n) j=1 if ((xmin * xmin + xmax * xmax ).lt. .0001) then xmin=y(1) xmax=y(n) endif do 200 i=1,n if((y(i).ge.xmin).and.(y(i).le.xmax)) then dy(j)=y(i) j=j+1 endif 200 continue n=j-1 write (*,*) 'Dataset :', dataset write (*,*) 'shifts : ', m write (*,*) 'Boundaries:', xmin,xmax write (*,*) 'Sample Size', n if(lambda.eq.0) lambda=log(dble(n)) / log(2.d0) write(*,*) 'Finest Resolution', lambda c xmax1=xmax delta=xmax-xmin dt=delta/dble(nplot - 1)*1.2d0 t(1)=xmin-delta*.1d0 fhad(1)=0.d0 do 300 i=2,nplot t(i)=t(i-1)+dt fhad(i)=0.d0 300 continue c do 800 i=1,m call bindenbu(dy,n,s,w,lambda,bb,nnb,xmin,xmax,nmax) call hist(dy,n,to,ny,nnb,bw,bb,nplot,t,fhad2) xmax=xmax+ delta/dble(m) do 900 j=1,nplot - 1 fhad(j) = fhad(j) + fhad2(j) 900 continue 800 continue xmax=xmax1 do 1000 i=1,m call bindenbu(dy,n,s,w,lambda,bb,nnb,xmin,xmax,nmax) call hist(dy,n,to,ny,nnb,bw,bb,nplot,t,fhad2) xmin=xmin-delta/dble(m) do 1100 j=1,nplot-1 fhad(j) = fhad(j) + fhad2(j) 1100 continue 1000 continue c do 400 i=1,nplot-1 write(50,500) t(i),fhad(i)/dble(2 * m) 400 continue fhad(nplot)=0.d0 close(50) close(60) close(70) 500 format(2f10.5) 600 format(f8.3) 700 format(2f8.3) c stop end subroutine parm(dataset,m, xmin,xmax,lambda) c c IMPLICIT undefined (a-z) real*8 xmin,xmax integer iargc, number,m,lambda character*30 arg, dataset c-----Ermittle die Anzahl der Parameter number=iargc() if ((number.gt.5).or.(number.lt.2)) then write (*,*) 'Subroutine parma: Wrong number of arguments.' write (*,*) 'Subroutine parma : Stop.' stop endif c-----Lese nacheinander die einzelnen Argumente call getarg(1, dataset) call getarg(2, arg) read (arg, *) m xmin=0.d0 xmax=0.d0 lambda=0 if (number.eq.5) then call getarg(3, arg) read (arg, *) xmin call getarg(4, arg) read (arg, *) xmax call getarg(5, arg) read (arg, *) lambda endif c c -----Fertig! return end subroutine bindenbu(y,n,s,w,lambda,b,nnb,xmin,xmax,nmax) cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c version February 1996 c c subroutine to calculate bins for rps density estimator c based on adaptive rule, bottom-up algorithm c c who is who c c y(n) input observations c n input number of observations (<=10000) c s (nmax) work score c w (nmax) work cell count c lambda input log_2 (n) (finest resolution level) c b(0:nmax) output bins c nnb output number of bins c xmin input left boundary c xmax input right boundary c nmax maximal size of data vector c c for details see: c J. Engel, c Adaptive Histograms based on Haar series c c subroutines: bsort,equalout c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc implicit undefined (a-z) integer i,j,k,lambda,m1,nb,n,nnb,nmax,w(nmax) double precision b(0:nmax),s(nmax),y(n),lhs,rhs,xn,xn1 double precision dx,xmin,xmax c xn=dble(n) xn1=2*xn/(xn+1) nb=2**lambda b(0)=xmin dx=xmax-xmin do 100 k=1,nb b(k)=xmin + dx* dble(k)/dble(2**lambda) !maximal partition w(k)=0 s(k)=0.d0 100 continue k=1 do 200 i=1,n 300 if(y(i).gt.b(k)) then k=k+1 goto 300 endif w(k)=w(k)+1 !cell counts 200 continue c m0=2**lambda do 400 j=1,lambda m1=2**(j-1) do 500 i=1,nb/m1,2 lhs=(w(i)-w(i+1))*(w(i)-w(i+1))+s(i) rhs= xn1*(w(i) + w(i+1)) w((i+1)/2)=w(i)+w(i+1) if (lhs.le.rhs) then do 600 k= m1*(i-1)+1,m1*(i+1)-1 b( k )=0.d0 600 continue c print*, 'knocking off b', 2**(j-1)*i else s((i+1)/2)=lhs-rhs endif c write(*,*) b(i),dmax(lhs-rhs,0) 500 continue 400 continue c call bsort(b,nb) call equalout(b,nb,b,nnb) c do i=1,nnb c write(*,*) b(i) c enddo c return end subroutine bsort(f,n) C------- SHELL SORT ALGORITHM CACM JULY 1964 implicit undefined (a-z) double precision f(n),x integer j,m,n,k,l m=n 20 m=m/2 k=n-m do 40 j=1,k l=j 50 if(f(l+m).ge.f(l)) goto 40 x=F(l+m) f(l+m)=f(l) f(l)=x if(l.le.m) goto 40 l=l-m goto 50 40 continue if(m.gt.1) goto 20 return end subroutine equalout(x,m,xx,n) c----------------------------------------- c c program to eliminate identical entries in c sorted data vector x c c x(0:m) input c m size of x c xx(0:n) output c n size of xx (<= m) c c------------------------------------------- implicit undefined (a-z) double precision x(0:m),xx(0:m) integer i,m, n n=0 xx(0)=x(0) do 1 i=1,m if (x(i).eq.xx(n)) goto 1 n=n+1 xx(n)=x(i) 1 continue return end subroutine hist(x,n,b0,ny,nnb,b,bb,m,t,f) C----------------------------------------------------------------------- C C VERSION: May 9 1993 C C PURPOSE: C C CALCULATES HISTOGRAM ESTIMATES C C PARAMETERS: C C INPUT X(N) DATA (Must BE SORTED) C INPUT N LENGTH OF X C INPUT B BIN SIZE C INPUT b0 ORIGIN(BETWEEN SMALLEST AND LARGEST DATA POINT) C INPUT nnb NUMBER OF BINS C INPUT M NUMBER OF POINTS TO ESTIMATE C INPUT NY 0: EQUAL CELL WIDTH C 1: VARIABLE CELL WIDTH C INPUT bb(K+1) SEQUENCE OF BIN BOUNDARY POINTS (FOR NY=1) C WORK bb(K+1) SEQUENCE OF BIN BOUNDARY POINTS (FOR NY=0) c INPUT b bin width (for ny=0) c C INPUT t(m) OUTPUT GRID (EQUIDISTANT) C C OUTPUT f(m) ESTIMATED DENSITY C C----------------------------------------------------------------------- implicit undefined (a-z) integer k,m,n,nnb,ny,i,j,jc double precision b0,b double precision x(n),f(m),bb(0:5*n),t(m) c c -- DEFINE BIN BOUNDARY POINTS IN CASE OF NY=0 c if (ny.eq.0) then bb(1)=b0 do 10 i=2,nnb bb(i)=bb(i-1)+b 10 continue endif c j=1 k=1 jc=0 i=1 c 100 continue if(x(i).le.bb(k)) then jc=jc+1 i=i+1 if(i.le.n) goto 100 endif 50 if((j.le.m).and.(t(j).le.bb(k))) then f(j)=jc/((bb(k)-bb(k-1))*n) j=j+1 goto 50 endif k=k+1 jc=0 if(i.le.n) goto 100 return end