c To use, execute commands: c c f77 -c modetree.f c ld -r -d -o mt.o modetree.o -lF77 -lm -lc c c then dyn.load mt.o into S c c April 8, 1986 c c Find bin counts of data array "x(n)" for ASH estimator c c "nbin" bins are formed over the interval [a,b) c c bin counts returned in array "bc" - # pts outside [a,b) = "nskip" c ##### Copyright 1986 David W. Scott subroutine bin1 ( x , n , ab , nbin , bc , nskip ) double precision x(n), bc(nbin), ab(2), a, b, d integer n, nbin, nskip, i, k nskip = 0 a = ab(1) b = ab(2) c call intpr("Into bin1",9,nskip,1) do 5 i = 1,nbin bc(i) = 0 5 continue d = (b-a) / nbin do 10 i = 1,n k = (x(i)-a) / d + 1.0 if (k.ge.1 .and. k.le.nbin) then bc(k) = bc(k) + 1 else nskip = nskip + 1 end if 10 continue return end c April 8, 1986 c c Computer ASH density estimate; Gaussian kernel c c Standard deviation of h c c Bin counts in array "bc(nbin)" - from routine "nbin1" c c "nbin" bins are formed over the interval [a,b) c c ASH estimates returned in array "f(nbin)" c c FP-ASH plotted at a+d/2 ... b-d/2 where d = (b-a)/nbin c c Note: If "nskip" was nonzero, ASH estimates incorrect near boundary c Note: Should leave "m" empty bins on each end of array "nc" so f OK c ##### Copyright 1986 David W. Scott c Modified to normal kernel 1991. Michael C. Minnotte subroutine nash1 ( h , bc , nbin , ab , t , f , w) double precision bc(nbin), ab(2), t(nbin), f(nbin), w(nbin), h, + pi, a, b, delta, c1, c2, cons, c integer nbin, ier, n, i, k pi = 3.1415926536 ier = 0 a = ab(1) b = ab(2) n = 0 delta = (b-a)/nbin c call intpr("Into nash1",10,ier,1) c-compute weights c w-array shifted by 1 c1 = 1/(h*sqrt(2*pi)) c2 = -1* delta**2/(2*h**2) w(1) = c1 cons = w(1) do 5 i = 2,nbin w(i) = c1*exp(c2*(i-1)**2) cons = cons + 2*w(i) 5 continue c do 6 i = 1, nbin c w(i) = w(i) / cons c6 continue c-compute ash(m) estimate do 10 i = 1,nbin t(i) = a + (i-0.5)*delta f(i) = 0.0 n = n + bc(i) 10 continue do 20 i = 1,nbin if (bc(i).eq.0) goto 20 c = bc(i) / n do 15 k = 1, nbin f(k) = f(k) + c * w( iabs(k-i)+1 ) 15 continue 20 continue return end c July 23, 1991 c c Finds the closest member of blist to a and returns its index as m1. c Finds closest member of blist on other side of a (if such exists) c and returns its index as m2. c nb is the number of elements of blist actually used. c ##### Copyright 1991 Michael C. Minnotte subroutine close(a, blist, nb, m1, m2) double precision a, blist(100), c, w integer nb, m1, m2, i c call intpr("Into close",10,nb,1) c = abs(a - blist(1)) + 1. do 10 i = 1, nb w = abs(a - blist(i)) if (w .lt. c) then m1 = i c = w end if 10 continue if ((a - blist(m1)) .gt. 0.) then m2 = m1 + 1 else if (m1 .eq. 1) then m2 = nb + 1 else m2 = m1 - 1 end if return end c July 23, 1991 c c Matches members of alist and blist (each array of length "len", filled to c "na" and "nb" respectively). Matched indices are output in mout(len,2). c "matwk1" and "matwk2" are work arrays (also (len,2)). c ##### Copyright 1991 Michael C. Minnotte subroutine matcher (alist, blist, na, nb, mout) double precision alist(100), blist(100) integer matwk1(100,2), matwk2(100,2), mout(100,2), na, nb, na1, + nb1, i c call intpr("Into matcher",12,na,1) na1 = na + 1 nb1 = nb + 1 do 10 i = 1, na mout(i,1) = nb1 call close (alist(i), blist, nb, matwk1(i,1), + matwk1(i,2)) 10 continue do 20 i = 1, nb mout(i,2) = na1 call close (blist(i), alist, na, matwk2(i,1), + matwk2(i,2)) 20 continue matwk1(na1,1) = nb1 matwk1(na1,2) = nb1 matwk2(nb1,1) = na1 matwk2(nb1,2) = na1 do 30 i = 1, na if (matwk2(matwk1(i,1),1) .eq. i) then mout(i,1) = matwk1(i,1) mout(matwk1(i,1),2) = i matwk2(matwk1(i,1),1) = na1 matwk2(matwk1(i,1),2) = na1 matwk1(i,1) = nb1 matwk1(i,2) = nb1 end if 30 continue do 40 i = 1, na if (matwk2(matwk1(i,1),2) .eq. i) then mout(i,1) = matwk1(i,1) mout(matwk1(i,1),2) = i matwk2(matwk1(i,1),1) = na1 matwk2(matwk1(i,1),2) = na1 matwk1(i,1) = nb1 matwk1(i,2) = nb1 end if 40 continue do 50 i = 1, na if (matwk2(matwk1(i,2),1) .eq. i) then mout(i,1) = matwk1(i,2) mout(matwk1(i,2),2) = i matwk2(matwk1(i,2),1) = na1 matwk2(matwk1(i,2),2) = na1 matwk1(i,1) = nb1 matwk1(i,2) = nb1 end if 50 continue do 60 i = 1, na if (matwk2(matwk1(i,2),2) .eq. i) then mout(i,1) = matwk1(i,2) mout(matwk1(i,2),2) = i matwk2(matwk1(i,2),1) = na1 matwk2(matwk1(i,2),2) = na1 matwk1(i,1) = nb1 matwk1(i,2) = nb1 end if 60 continue return end c July 23, 1991 c c Finds the vertex of the parabola going through the points c (x1, y1), (x2, y2), and (x3, y3). c c ##### Copyright 1991 Michael C. Minnotte double precision function quad (x1, x2, x3, y1, y2, y3) double precision x1, x2, x3, y1, y2, y3, a, b c call intpr("Into quad",9,0,1) a = ((y1-y2)/(x1-x2) - (y2-y3)/(x2-x3))/(x1-x3) b = (y1-y2)/(x1-x2) - a*(x1+x2) quad = -b/(2.*a) return end c July 23, 1991 c c Finds the modes of a function f(nbin) on t(nbin). Finds each c triplet of values in which the middle is highest of the three, c then uses quad (above) to estimate the actual mode. Returns a c list of modes (rm), and the indices of the antimodes on either side c of each mode (ml and mr), along with the total number of modes (nm). c Currently limited; will not find modes with plateaus instead of sharp c peaks. May need to be improved to do this. c c ##### Copyright 1991 Michael C. Minnotte subroutine findmodes (f, t, nbin, rm, mc, ml, mr, nm) double precision f(nbin), t(nbin), rm(100), quad integer mc(100), ml(100), mr(100), nbin, nm, i, j nm = 0 c call intpr("Into findmodes",14,nm,1) do 5 i = 1, 100 mr(i) = 0 5 continue ml(1) = 1 do 10 i = 2, (nbin-1) if ((f(i) .gt. f(i-1)) .and. (f(i) .ge. f(i+1))) then nm = nm + 1 mc(nm) = i rm(nm) = quad(t(i-1),t(i),t(i+1),f(i-1),f(i),f(i+1)) else if ((f(i) .lt. f(i-1)) .and. (f(i) .le. f(i+1))) then mr(nm) = i ml(nm + 1) = i else if ((f(i) .le. f(i-1)) .and. (f(i) .lt. f(i+1))) then ml(nm + 1) = i end if 10 continue mr(nm) = nbin i = 1 20 if (i .gt. (nm-1)) goto 40 if (mr(i).eq.0) then mc(i) = mc(i+1) mr(i) = mr(i+1) rm(i) = rm(i+1) do 30 j = (i+1), (nm-1) ml(j) = ml(j+1) mc(j) = mc(j+1) mr(j) = mr(j+1) rm(j) = rm(j+1) 30 continue nm = nm - 1 else i = i+1 end if goto 20 40 continue return end c July 28, 1993 c c takes density f(nbin) and returns it in g(nbin) with points with c indices between ml and mr replaced with the smaller of f(ml) and f(mr) c d (passed) is the bin widths; coeff (returned) is the total area c under the curve for g. subroutine excise (f, g, coeff, nbin, ml, mr, d) double precision f(nbin), g(nbin), coeff, d, tot, val integer nbin, ml, mr, incrmt, i, j tot = 0.d0 if (f(ml) .gt. f(mr)) then incrmt = 1 val = f(ml) i = ml + 1 do 10 j = 1, ml g(j) = f(j) tot = tot + f(j) 10 continue else incrmt = -1 val = f(mr) i = mr - 1 do 20 j = mr, nbin g(j) = f(j) tot = tot + f(j) 20 continue end if 30 if (f(i) .le. val) goto 40 g(i) = val tot = tot + val i = i + incrmt goto 30 40 continue 50 if ((i .eq. 0) .or. (i .eq. (nbin+1))) goto 60 g(i) = f(i) tot = tot + f(i) i = i + incrmt goto 50 60 continue coeff = tot * d c Uncomment the following lines to rescale the new density properly c (Leave commented for enhanced mode tree) c do 70 i = 1, nbin c g(i) = g(i)/coeff c 70 continue return end