program mrshis cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Version from February 3, 1996 c programmer: Joachim Engel c c program for multiresolution histogram c c Call: mrshis dataset xmin xmax lambda c c dataset set of univariate data (may be unordered) c xmin c xmax algorithm applies to interval [xmin,xmax] c lambda finest resolution level c c e.g. mrshis faithful 1.2 5. 5 c default for boundaries: smallest and largest observation c resolution level: log(sample size) c c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc implicit undefined (a-z) c integer nmax,n,nnb,nplot,ny,i,j,lambda character*30 dataset parameter(nplot=2000,nmax=10000) integer w(nmax) double precision t(nplot),fhad(nplot) double precision dy(nmax),y(nmax),bb(nmax),s(nmax) double precision xmin,xmax,dt,bw,to data ny/1/,to/0/ c call parm(dataset,xmin,xmax,lambda) c open(unit=50,file='mrshis.est') open(unit=60,file=dataset) open(unit=70,file='tmp.dat') do 100 i=1,nmax read(60,*,end=3) y(i) write(70,700) y(i),0.d0 100 continue 3 n=i-1 call bsort(y,n) c do i=1,n c print*, i, y(i) c enddo 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 (*,*) 'Boundaries:', xmin,xmax write (*,*) 'Sample Size', n c dt=(xmax-xmin)/dble(nplot) t(1)=xmin do 300 i=2,nplot t(i)=t(i-1)+dt 300 continue if (lambda.eq.0) lambda=log(dble(n)) / log(2.d0) -1 c lambda=log(dble(n)) / log(2.d0) -1 write (*,*) 'Finest Resolution ', lambda c call bindenbu(dy,n,s,w,lambda,bb,nnb,xmin,xmax,nmax) call hist(dy,n,to,ny,nnb,bw,bb,nplot,t,fhad) c do 400 i=1,nplot write(50,500) t(i),fhad(i) 400 continue close(50) close(60) close(70) 500 format(2f10.5) 600 format(f8.3) 700 format(2f8.3) c stop end subroutine parm(dataset, xmin,xmax,lambda) c c IMPLICIT undefined (a-z) real*8 xmin,xmax integer iargc, number,lambda character*30 arg, dataset c-----Ermittle die Anzahl der Parameter number=iargc() if ((number.gt.4).or.(number.lt.1)) then write (*,*) 'Subroutine parma: Wrong number of arguments.' write (*,*) 'Subroutine parma : Stop.' stop endif c-----Lese nacheinander die einzelnen Argumente call getarg(1, dataset) xmin=0.d0 xmax=0.d0 lambda=0 if (number.eq.4) then call getarg(2, arg) read (arg, *) xmin call getarg(3, arg) read (arg, *) xmax call getarg(4, arg) read (arg, *) lambda endif c c -----Fertig! 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 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 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:3*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