cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c c
c Fortran Mode Testing Software c
c c
c Save as text under tme.f, then compile as: c
c f77 -u -C -o tme tme.f /usr/local/lib/libnag.a c
c c
c If NAG library is unavailable, replace G05* calls below with another c
c random number generator. c
c c
c Calculates mode tree, performs individual mode tests as in c
c Minnotte, M.C., (1997), "Nonparametric testing of the c
c existence of modes," The Annals of Statistics, 25, 1646-1660. c
c c
c Copyright 1995, Michael C. Minnotte c
c c
c March 28, 1995 c
c c
c Data should be one-to-a-line in file 'dataxyz'. c
c Output will be found in files 'treeoutxyz' and 'modeoutxyz'. c
c c
c See associated S-plus code for easy input/output, and plotting of c
c resulting tree. c
c c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
program tme
integer i, j, k, n, nbin, maxx, nskip, nh, count,
+ mco(100), mcn(100), mlo(100), mln(100),
+ mro(100), mrn(100), nmo, nmn, mout(100,2),
+ oldfrom, newfrom, npass(2), nbig, nmax, nm
parameter (nbin=500, maxx=5000, nh=200, nbig=16, nmax=400)
double precision x(maxx), bins(nbin),
+ ab(2), delta, hlist(nh), dh, t(nbin),
+ f(nbin), w(nh,nbin), rmo(100), rmn(100), pval,
+ bbc(nbin), g(nbin), fb(nbin), gb(nbin), wk(nbin),
+ probs(nbin), rmode, scale, cp(nbin)
parameter (scale = 1D100)
open(unit=2,file='dataxyz')
open(unit=3,file='treeoutxyz')
open(unit=4,file='modeoutxyz')
c read in data, find min, max
n = 1
read(2,*,end=120) x(1)
ab(1) = x(1)
ab(2) = x(2)
do 110 i = 2, maxx
read(2,*,end=120) x(i)
if (x(i).lt.ab(1)) then
ab(1) = x(i)
else if (x(i).gt.ab(2)) then
ab(2) = x(i)
end if
c write (*,*) 'Point ',i,' = ',x(i)
n = n + 1
110 continue
120 continue
c write (*,*) 'n = ',n
c extend range of data 10% each direction
delta = (ab(2)-ab(1))*0.1
ab(1) = ab(1) - delta
ab(2) = ab(2) + delta
c bin data
call bin(x, maxx, n, ab, nbin, bins, nskip)
c choose bandwidths
hlist(1) = ((ab(2)-ab(1))/4)
hlist(nh) = ((ab(2)-ab(1))/200)
if (nh .eq. 1) then
write (*,*) 'nh = 1. Div by 0 in main.'
end if
delta = (log(hlist(1)) - log(hlist(nh)))/(nh-1)
do 210 i = 2, (nh-1)
hlist(i) = exp(log(hlist(1)) - (i-1)*delta)
210 continue
if (hlist(1) .lt. 1D-50) then
write (*,*) 'hlist(1) = ',hlist(1),'. Div by 0 in main.'
end if
dh = hlist(2)/hlist(1)
npass(1) = nbig
npass(2) = nmax
c first estimate
call centers(t, ab, nbin, delta)
c G05CBF(0) - repeatable initialize of random number generator (NAG lib.)
c G05CCF - nonrepeatable initialize
c if NAG unavailable, replace with initialization (if any) of rng of choice
call G05CBF(0)
call kweights(hlist, w, nbin, nh, n, delta, scale)
call nash(bins, nbin, f, w, nh, 1)
call findmodes(f, t, nbin, rmo, mco, mlo, mro, nmo)
c first tree records
do 310 i = 1, nmo
write (3,3000) rmo(i), hlist(1), 0, 0
310 continue
count = 0
c loop over h
do 1010 i = 2, nh
write (*,*) 'h number', i, ' starting.'
call nash(bins, nbin, f, w, nh, i)
call findmodes(f, t, nbin, rmn, mcn, mln, mrn, nmn)
call matcher(rmn, rmo, nmn, nmo, mout)
c loop over modes
do 510 j = 1, nmn
if (mout(j,1).le.nmo) then
write (3,3000) rmn(j), hlist(i), count + mout(j,1), 1
else
do 410 k = 1, nmo
if (mlo(k).lt.mcn(j) .and. mro(k).gt.mcn(j)) then
oldfrom = k
end if
410 continue
newfrom = mout(oldfrom, 2)
write (3,3000) rmn(j), hlist(i), count + nmo + newfrom,
+ 2
if (nmo.eq.1) then
write (4,4000) rmo(1), hlist(i-1), 0, t(mlo(1)),
+ t(mro(1))
else
call modetest(bins, bbc, f, g, fb, gb, t, wk,
+ probs, nbin, rmode, pval, hlist,
+ oldfrom, nm, npass, n, nh, i, delta, w, cp,
+ scale)
write (4,4000) rmo(oldfrom), hlist(i-1), pval,
+ t(mlo(oldfrom)), t(mro(oldfrom))
write (*,*) rmo(oldfrom), hlist(i-1), pval,
+ t(mlo(oldfrom)), t(mro(oldfrom))
end if
end if
510 continue
count = count + nmo
nmo = nmn
do 610 j = 1, nmn
rmo(j) = rmn(j)
mco(j) = mcn(j)
mlo(j) = mln(j)
mro(j) = mrn(j)
610 continue
1010 continue
3000 format (F12.6, F12.6, I6, I4)
4000 format (5F12.6)
c i = ieee_flags('clear','exception','all')
stop
end
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 bin ( x , maxx , n , ab , nbin , bc , nskip )
double precision x(maxx), bc(nbin), ab(2), a, b, d
integer n, nbin, nskip, i, k, maxx
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
if (nbin .eq. 0) then
write (*,*) 'nbin = 0. Div by 0 in bin.'
end if
d = (b-a) / nbin
if (d .lt. 1D-50) then
write (*,*) 'd = ',d,'. Div by 0 in bin.'
end if
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 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 24, 1991
c
c The biggie! Finds the modes of the ash of bc with nbin bins and
c smoothing over m bins. For each mode gets a p-value by removing it,
c then finding the number of cases, out of nb bootstraps from the
c modified distribution, in which there is a matching mode, at least as
c "extreme" (defined for now as ise between the two defining antimodes)
c as the one originally excised. Key output is nm (number of modes),
c rmodes (first nm elements are the modes), and pval (first nm elements
c are corresponding p-values). Input includes bc, nbin, dh,
c ab, and m. All else is work arrays for various subroutines.
c
c Copyright 1991 Michael C. Minnotte
subroutine modetest(bc, bbc, f, g, fb, gb, t, wk, probs, nbin,
+ rmode, pval, hlist, k, nm, npass, n, nh, ih, delta, w,
+ cp, scale)
double precision bc(nbin), bbc(nbin), f(nbin), g(nbin),
+ gb(nbin), t(nbin), wk(nbin), probs(nbin), rmodes(100),
+ rmode, pval, hlist(nh), delta, coeff, exr, w(nh, nbin),
+ exb, fb(nbin), aid, cp(nbin), scale
integer nbin, k, nm, nbig, nmax, ip1, ip2, nmax1, n, ncount,
+ ml(100), mr(100), npass(2), j, mc(100), nh, ih
c call intpr("Into modetest",13,k,1)
nbig = npass(1)
nmax = npass(2)
call nash(bc, nbin, f, w, nh, ih)
call findmodes(f, t, nbin, rmodes, mc, ml, mr, nm)
c call inflpt(f, mc(k-1), mr(k-1), nbin, ip1)
c call inflpt(f, ml(k+1), mc(k+1), nbin, ip2)
rmode = rmodes(k)
nmax1 = nmax - 1
call excise(f, g, coeff, nbin, ml(k), mr(k), delta)
pval = 0.0
exr = aid(f, g, ml(k), mr(k), nbin)
c call dblepr("exr",3,exr,1)
call restore(f,g,wk,nbin,delta,mc,ml,mr,nm,k,scale)
ncount = 0
j = 1
10 continue
c if (j .eq. 1) call intpr("simmode call #",14,j,1)
c edit to pass correct values to simmode, return exb
call simmode(hlist,k,g,rmodes,fb,bbc,t,wk,nbin,delta,
+ n,nm,ip1,ip2,exb,ih,w,nh,cp,probs,scale)
if (exb .ge. exr) ncount = ncount + 1
if (ncount .eq. nbig) then
if (j .eq. 0) then
write (*,*) 'j = 0. Div by 0 in modetest.'
end if
pval = real(nbig)/real(j)
else if (j .eq. nmax1) then
if (nmax .eq. 0) then
write (*,*) 'nmax = 0. Div by 0 in modetest.'
end if
pval = real(ncount+1)/real(nmax)
else
j = j + 1
goto 10
end if
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
if (nm .ge. 1) then
mr(nm) = i
ml(nm + 1) = i
end if
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 ***************************************************************
subroutine kweights (hlist, w, nbin, nh, n, delta, scale)
integer nbin, nh, n, i, j
double precision hlist(nh), w(nh,nbin), c1, c2, delta,
+ scale, pi
pi = 3.1415926536
do 20 j = 1, nh
if (hlist(j)*sqrt(2*pi)*n .lt. 1D-50) then
write (*,*) 'Div by 0 in kweights. c1.'
end if
if (2*hlist(j)**2 .lt. 1D-50) then
write (*,*) 'Div by 0 in kweights. c2.'
end if
c1 = scale/(hlist(j)*sqrt(2*pi)*n)
c2 = -1* delta**2/(2*hlist(j)**2)
w(j,1) = c1
do 10 i = 2,nbin
w(j,i) = c1*exp(c2*(i-1)**2)
10 continue
20 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 nash ( bc , nbin , f , w, nh, ih)
double precision bc(nbin), f(nbin), w(nh,nbin)
integer nbin, i, k, nh, ih
c-compute ash(m) estimate
do 10 i = 1,nbin
f(i) = 0.0
10 continue
do 20 i = 1,nbin
if (bc(i).gt.0) then
do 15 k = 1, nbin
f(k) = f(k) + bc(i) * w(ih, iabs(k-i)+1 )
15 continue
end if
20 continue
return
end
c ******************************************************************
subroutine centers(t, ab, nbin, delta)
integer nbin, i
double precision ab(2), t(nbin), delta
if (nbin .eq. 0) then
write(*,*) 'nbin = 0. Div by 0 in centers.'
end if
delta = (ab(2) - ab(1))/nbin
do 10 i = 1, nbin
t(i) = ab(1) + (i-0.5)*delta
10 continue
return
end
c July 23, 1991
c
c Provides an (unscaled) estimate of absolute integrated distance between
c functions f and g, recorded over the same nbin equally spaced points.
c Calculates AID over bins ml to mr (inclusive).
c
c ##### Copyright 1991 Michael C. Minnotte
double precision function aid (f, g, ml, mr, nbin)
double precision f(nbin), g(nbin), raid
integer ml, mr, nbin, i
c call intpr("Into aid",8,ml,1)
raid = 0.d0
do 10 i = ml, mr
raid = raid + abs(f(i) - g(i))
10 continue
aid = raid
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 April 21, 1992
c simulate from g,
subroutine simmode(hlist, k, g, rmodes, fb, bbc, t, wk,
+ nbin, delta, n, nm, ip1, ip2, exb, ih, w, nh,
+ cp, probs, scale)
double precision g(nbin), fb(nbin), bbc(nbin), t(nbin),
+ wk(nbin), probs(nbin), hlist(nh), delta, exb, rmdmax, cb,
+ rmodes(100), bmodes1(100), bmodes2(100), aid,
+ w(nh,nbin), cp(nbin), scale
integer mcb(100), mlb(100), mrb(100), mout(100,2), k, nbin, n,
+ nm, ip1, ip2, mli, mri, mdmax, i, nmb1, nmb2, nwmmx, nh, ih,
+ ih2
c call intpr("Into simmode",12,k,1)
call cprobs(g, probs, cp, nbin, delta, scale)
call sim(bbc, cp, nbin, n, scale)
call nash(bbc, nbin, fb, w, nh, ih)
call findmodes(fb, t, nbin, bmodes1, mcb, mlb, mrb, nmb1)
call matcher(rmodes, bmodes1, nm, nmb1, mout)
c find range where match may be found
if (k .eq. 1) then
mli = 1
else if (mout(k-1,1) .le. nmb1) then
mli = mrb(mout(k-1,1)) + 1
else
mli = ip1 + 1
end if
if (k .eq. nm) then
mri = nbin
else if (mout(k+1,1) .le. nmb1) then
mri = mlb(mout(k+1,1)) - 1
else
mri = ip2 - 1
end if
c call intpr("mli",3,mli,1)
c call intpr("mri",3,mri,1)
c find largest mode (if any) within range
rmdmax = 0.0
mdmax = 0
do 10 i = 1, nmb1
if (mcb(i) .lt. mli) goto 10
if (mcb(i) .gt. mri) goto 10
call excise(fb, wk, cb, nbin, mlb(i), mrb(i), delta)
exb = aid(fb, wk, mlb(i), mrb(i), nbin)
if (exb .gt. rmdmax) then
mdmax = i
rmdmax = exb
end if
10 continue
c skips decreasing h section; remove ?????
c call intpr("mdmax",5,mdmax,1)
c call dblepr("rmdmax",6,rmdmax,1)
c exb = rmdmax
c goto 500
c end of optional part
c find h for which mode is largest
if (mdmax .eq. 0) then
exb = 0
goto 500
else if (ih .eq. nh) then
exb = rmdmax
goto 500
else
ih2 = ih
20 ih2 = ih2 + 1
call nash(bbc, nbin, fb, w, nh, ih2)
call findmodes(fb, t, nbin, bmodes2, mcb, mlb, mrb, nmb2)
call matcher(bmodes1, bmodes2, nmb1, nmb2, mout)
nwmmx = mout(mdmax,1)
c call intpr("nwmmx",5,nwmmx,1)
c call intpr("mlb(nwmmx)",10,mlb(nwmmx),1)
c call intpr("mrb(nwmmx)",10,mrb(nwmmx),1)
if (nwmmx .gt. nmb2) then
exb = rmdmax
c call dblepr("h2",2,h2,1)
c call intpr("ic",2,ic,1)
c call dblepr("exb",3,exb,1)
goto 500
end if
call excise(fb, wk, cb, nbin, mlb(nwmmx), mrb(nwmmx), delta)
exb = aid(fb, wk, mlb(nwmmx), mrb(nwmmx), nbin)
if ((exb .gt. rmdmax) .and. (ih2 .lt. nh)) then
rmdmax = exb
mdmax = nwmmx
nmb1 = nmb2
do 30 i = 1, 100
bmodes1(i) = bmodes2(i)
30 continue
goto 20
else
exb = rmdmax
c call dblepr("h2",2,h2,1)
c call intpr("ic",2,ic,1)
c call dblepr("exb",3,exb,1)
goto 500
end if
end if
500 return
end
c April 16, 1992
c find L_1 approximation to f without mode k
c
subroutine restore(f,g,fr,nbin,d,mc,ml,mr,nm,k,scale)
double precision f(nbin), g(nbin), fr(nbin), d, scale
integer mc(100), ml(100), mr(100), mcr(100), mlr(100), mrr(100),
+ k, nbin, nm, i
c call intpr("Into restore",12,k,1)
c call intpr("nm",2,nm,1)
if (k .eq. 1) then
do 10 i = 1, nbin
g(i) = f(i)
10 continue
c call intpr("Calling rest_end",16,1,1)
call rest_end(f,g,nbin,d,mc,ml,mr,scale)
else if (k .eq. nm) then
c call intpr("Calling reverse",15,nm,1)
call reverse(f,g,nbin)
do 20 i = 1, nm
mcr(i) = nbin + 1 - mc(nm + 1 - i)
mlr(i) = nbin + 1 - mr(nm + 1 - i)
mrr(i) = nbin + 1 - ml(nm + 1 - i)
20 continue
do 30 i = 1, nbin
fr(i) = g(i)
30 continue
c call intpr("Calling rest_end",16,nm,1)
call rest_end(g,fr,nbin,d,mcr,mlr,mrr,scale)
call reverse(fr,g,nbin)
else
call rest_mid(f,g,nbin,d,mc,ml,mr,k,scale)
end if
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)
if (abs(x1-x2) .lt. 1D-50) then
write (*,*) 'x1-x2 = ',x1-x2,'. Div by 0 in quad.'
end if
if (abs(x2-x3) .lt. 1D-50) then
write (*,*) 'x2-x3 = ',x2-x3,'. Div by 0 in quad.'
end if
if (abs(x1-x3) .lt. 1D-50) then
write (*,*) 'x1-x3 = ',x1-x3,'. Div by 0 in quad.'
end if
a = ((y1-y2)/(x1-x2) - (y2-y3)/(x2-x3))/(x1-x3)
b = (y1-y2)/(x1-x2) - a*(x1+x2)
if (abs(a) .lt. 1D-50) then
write (*,*) 'a = ',a,'. Div by 0 in quad.'
end if
quad = -b/(2.*a)
return
end
subroutine excise (f, g, coeff, nbin, ml, mr, d)
double precision f(nbin), g(nbin), coeff, d, val, tot
integer nbin, ml, mr, incrmt, i, j
c call intpr("Into excise",11,ml,1)
tot = 0.d0
if (f(ml) .gt. f(mr)) then
incrmt = 1
val = f(ml)
i = ml + 1
c iflagl = 0
c iflagr = 0
c call intpr("In 1st loop",11,i,1)
c call intpr("ml",2,ml,1)
c call intpr("mr",2,mr,1)
do 10 j = 1, ml
g(j) = f(j)
tot = tot + f(j)
c if ((ml .gt. nbin) .and. (iflagl .eq. 0)) then
c call intpr("ml",2,ml,1)
c call intpr("j",1,j,1)
c iflagl = 1
c end if
c if ((mr .gt. nbin) .and. (iflagr .eq. 0)) then
c call intpr("mr",2,mr,1)
c call intpr("j",1,j,1)
c iflagr = 1
c end if
10 continue
else
incrmt = -1
val = f(mr)
i = mr - 1
c iflagl = 0
c iflagr = 0
c call intpr("In 2nd loop",11,i,1)
c call intpr("ml",2,ml,1)
c call intpr("mr",2,mr,1)
do 20 j = mr, nbin
g(j) = f(j)
tot = tot + f(j)
c if ((ml .gt. nbin) .and. (iflagl .eq. 0)) then
c call intpr("ml",2,ml,1)
c call intpr("j",1,j,1)
c iflagl = 1
c end if
c if ((mr .gt. nbin) .and. (iflagr .eq. 0)) then
c call intpr("mr",2,mr,1)
c call intpr("j",1,j,1)
c iflagr = 1
c end if
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 do 70 i = 1, nbin
c g(i) = g(i)/coeff
c 70 continue
return
end
c April 16, 1992
c reverse f, putting it into fr
c
subroutine reverse(f,fr,nbin)
double precision f(nbin), fr(nbin)
integer nbin, i
c call intpr("Into reverse",nbin,1)
do 10 i = 1, nbin
fr(i) = f(nbin + 1 - i)
10 continue
return
end
c ***************************************************************
subroutine cprobs(g, probs, cp, nbin, delta, scale)
integer nbin, i
double precision g(nbin), probs(nbin), cp(nbin), delta,
+ scale, tot
tot = 0.0
do 10 i = 1, nbin
probs(i) = g(i) * delta
tot = tot + probs(i)
10 continue
if (tot .lt. 1D-50) then
write (*,*) 'tot = ',tot,'. Div by 0 in cprobs.'
end if
tot = scale/tot
cp(1) = probs(1) * tot
do 20 i = 2, nbin
cp(i) = cp(i-1) + probs(i)*tot
20 continue
return
end
c July 23, 1991
c
c Simulates data from distribution g(nbin).
c Sample will return binned in bbc(nbin). d is the binwidth, and n is the
c total number of points to be sampled.
c
c ##### Copyright 1991 Michael C. Minnotte
subroutine sim (bbc, cp, nbin, n, scale)
double precision bbc(nbin), cp(nbin), u, scale,
+ G05CAF, dummy
integer nbin, n, i,ilow, ihigh, itest
c call intpr("Into sim",8,nbin,1)
do 10 i = 1, nbin
bbc(i) = 0.0
10 continue
c If NAG library unavailable, replace G05CAF(dummy) with another call to
c a uniform(0,1) random number generator.
do 30 i = 1, n
u = G05CAF(dummy)*scale
if (cp(1) .ge. u) then
bbc(1) = bbc(1) + 1
else
ilow = 1
ihigh = nbin
25 itest = (ihigh + ilow) / 2
if (cp(itest) .lt. u) then
ilow=itest
else
ihigh=itest
end if
if (ihigh-ilow .gt. 1) goto 25
bbc(ihigh) = bbc(ihigh) + 1
end if
30 continue
return
end
c April 16, 1992
c restore mass of left-most mode of density f
c
subroutine rest_end(f,g,nbin,d,mc,ml,mr,scale)
double precision f(nbin), g(nbin), d, rmax, rmin, rtest, area,
+ a1, a2,scale
integer mc(100), ml(100), mr(100), nbin, i
c call intpr("Into rest_end",13,nbin,1)
if (f(mc(2)) .lt. f(mc(1))) then
call checklev(f,g,nbin,mc,ml,mr,1,1,f(mc(2)),a1,a2)
if (a1 .ge. a2) goto 500
rmax = f(mc(2))
else
rmax = f(mc(1))
end if
rmin = f(mr(1))
call findcut(f,g,nbin,mc,ml,mr,1,1,rmax,rmin,rtest)
500 continue
area = 0
do 510 i = 1, nbin
area = area + g(i)
510 continue
if (scale .lt. 1D-50) then
write (*,*) 'scale = ',scale,'. Div by 0 in rest_end.'
end if
area = area*d/scale
if (area .lt. 1D-50) then
write (*,*) 'area = ',area,'. Div by 0 in rest_end.'
end if
do 520 i = 1, nbin
g(i) = g(i) / area
520 continue
return
end
c May 26, 1992
c restore mass of interior modes
c
subroutine rest_mid(f,g,nbin,d,mc,ml,mr,k,scale)
double precision f(nbin), g(nbin), d, rlev1, rlev2, a1, a2,
+ diff, area, scale
integer mc(100), ml(100), mr(100), nbin, k, mlout,
+ mrout, i
c call intpr("Into rest_mid",13,k,1)
if (f(mr(k)) .ge. f(mc(k-1))) then
if (f(mc(k)) .gt. f(mc(k+1))) then
call checklev(f,g,nbin,mc,ml,mr,k,1,f(mc(k+1)),a1,a2)
if (a1 .gt. a2) then
rlev2 = f(mc(k+1))
diff = a1 - a2
call checklev(f,g,nbin,mc,ml,mr,k,-1,f(mc(k-1)),a1,a2)
if (a2 .gt. diff) then
call findcut2(f,g,nbin,mc,ml,mr,k,-1,f(mc(k-1)),
+ f(ml(k)),diff,rlev1)
else
rlev1 = f(mc(k-1))
end if
call level(f, rlev1, nbin, mc(k-1), mc(k), ml(k),
+ mlout, mrout)
call checklev(f,g,nbin,mc,ml,mr,k,1,rlev2,a1,a2)
do 10 i = mlout, mrout
g(i) = rlev1
10 continue
else
call findcut(f,g,nbin,mc,ml,mr,k,1,f(mc(k+1)),f(mr(k)),
+ rlev2)
end if
else
call findcut(f,g,nbin,mc,ml,mr,k,1,f(mc(k)),f(mr(k)),
+ rlev2)
end if
else if (f(ml(k)) .ge. f(mc(k+1))) then
if (f(mc(k)) .gt. f(mc(k-1))) then
call checklev(f,g,nbin,mc,ml,mr,k,-1,f(mc(k-1)),a1,a2)
if (a1 .gt. a2) then
rlev1 = f(mc(k-1))
diff = a1 - a2
call checklev(f,g,nbin,mc,ml,mr,k,1,f(mc(k+1)),a1,a2)
if (a2 .gt. diff) then
call findcut2(f,g,nbin,mc,ml,mr,k,1,f(mc(k+1)),
+ f(mr(k)),diff,rlev2)
else
rlev2 = f(mc(k+1))
end if
call level(f, rlev2, nbin, mc(k), mc(k+1), mr(k),
+ mlout, mrout)
call checklev(f,g,nbin,mc,ml,mr,k,-1,rlev1,a1,a2)
do 20 i = mlout, mrout
g(i) = rlev2
20 continue
else
call findcut(f,g,nbin,mc,ml,mr,k,-1,f(mc(k-1)),
+ f(ml(k)),rlev1)
end if
else
call findcut(f,g,nbin,mc,ml,mr,k,-1,f(mc(k)),f(ml(k)),
+ rlev1)
end if
else
if (f(mr(k)) .gt. f(ml(k))) then
call checklev(f,g,nbin,mc,ml,mr,k,-1,f(mr(k)),a1,a2)
if (a1 .gt. a2) then
if (f(mc(k-1)) .gt. f(mc(k))) then
call findcut(f,g,nbin,mc,ml,mr,k,-1,f(mc(k)),
+ f(mr(k)),rlev1)
else
call checklev(f,g,nbin,mc,ml,mr,k,-1,f(mc(k-1)),
+ a1,a2)
if (a1 .gt. a2) then
rlev1 = f(mc(k-1))
else
call findcut(f,g,nbin,mc,ml,mr,k,-1,f(mc(k-1)),
+ f(mr(k)),rlev1)
end if
end if
else
call findcut2(f,g,nbin,mc,ml,mr,k,-1,f(mr(k)),
+ f(ml(k)),a1,rlev1)
end if
if (f(mc(k+1)) .gt. f(mc(k))) then
call findcut(f,g,nbin,mc,ml,mr,k,1,f(mc(k)),
+ f(mr(k)),rlev2)
else
call checklev(f,g,nbin,mc,ml,mr,k,1,f(mc(k+1)),
+ a1,a2)
if (a1 .gt. a2) then
rlev2 = f(mc(k+1))
else
call findcut(f,g,nbin,mc,ml,mr,k,1,f(mc(k+1)),
+ f(mr(k)),rlev2)
end if
end if
else
call checklev(f,g,nbin,mc,ml,mr,k,1,f(ml(k)),a1,a2)
if (a1 .gt. a2) then
if (f(mc(k+1)) .gt. f(mc(k))) then
call findcut(f,g,nbin,mc,ml,mr,k,1,f(mc(k)),
+ f(ml(k)),rlev2)
else
call checklev(f,g,nbin,mc,ml,mr,k,1,f(mc(k+1)),
+ a1,a2)
if (a1 .gt. a2) then
rlev2 = f(mc(k+1))
else
call findcut(f,g,nbin,mc,ml,mr,k,1,f(mc(k+1)),
+ f(ml(k)),rlev2)
end if
end if
else
call findcut2(f,g,nbin,mc,ml,mr,k,1,f(ml(k)),
+ f(mr(k)),a1,rlev2)
end if
if (f(mc(k-1)) .gt. f(mc(k))) then
call findcut(f,g,nbin,mc,ml,mr,k,-1,f(mc(k)),
+ f(ml(k)),rlev1)
else
call checklev(f,g,nbin,mc,ml,mr,k,-1,f(mc(k-1)),
+ a1,a2)
if (a1 .gt. a2) then
rlev1 = f(mc(k-1))
else
call findcut(f,g,nbin,mc,ml,mr,k,-1,f(mc(k-1)),
+ f(ml(k)),rlev1)
end if
end if
end if
if (rlev1 .gt. rlev2) then
call checklev(f,g,nbin,mc,ml,mr,k,-1,rlev1,a1,a2)
if (a1 .lt. 1D-50) then
write(*,*) 'a1 = ',a1,'. Div by 0 in rest_mid'
end if
if ((abs(a1-a2)/a1) .lt. 0.01) then
goto 500
else
diff = a1 - a2
call checklev(f,g,nbin,mc,ml,mr,k,1,rlev2,a1,a2)
if (a2 .gt. diff) then
call findcut2(f,g,nbin,mc,ml,mr,k,1,rlev2,f(mr(k)),
+ diff,rlev2)
end if
call level(f,rlev2,nbin,mc(k),mc(k+1),mr(k),mlout,
+ mrout)
call checklev(f,g,nbin,mc,ml,mr,k,-1,rlev1,a1,a2)
do 30 i = mlout, mrout
g(i) = rlev2
30 continue
goto 500
end if
else
call checklev(f,g,nbin,mc,ml,mr,k,1,rlev2,a1,a2)
if (a1 .lt. 1D-50) then
write(*,*) 'a1 = ',a1,'. Div by 0 in rest_mid'
end if
if ((abs(a1-a2)/a1) .lt. 0.01) then
goto 500
else
diff = a1 - a2
call checklev(f,g,nbin,mc,ml,mr,k,-1,rlev1,a1,a2)
if (a2 .gt. diff) then
call findcut2(f,g,nbin,mc,ml,mr,k,-1,rlev1,f(ml(k)),
+ diff,rlev1)
end if
call level(f,rlev1,nbin,mc(k-1),mc(k),ml(k),mlout,
+ mrout)
call checklev(f,g,nbin,mc,ml,mr,k,1,rlev2,a1,a2)
do 40 i = mlout, mrout
g(i) = rlev1
40 continue
goto 500
end if
end if
end if
500 continue
area = 0
do 510 i = 1, nbin
area = area + g(i)
510 continue
if (scale .lt. 1D-50) then
write (*,*) 'scale = ',scale,'. Div by 0 in rest_mid.'
end if
area = area*d/scale
if (area .lt. 1D-50) then
write (*,*) 'area = ',area,'. Div by 0 in rest_mid.'
end if
do 520 i = 1, nbin
g(i) = g(i) / area
520 continue
return
end
c April 16, 1992
c find level to get cut area = diff
subroutine findcut2(f,g,nbin,mc,ml,mr,k,idir,rmaxi,rmini,
+ diff, rtest)
double precision f(nbin), g(nbin), rmaxi, rmini, diff, rtest,
+ rmax, rmin, a1, a2
integer mc(100), ml(100), mr(100), nbin, k, idir
c call intpr("Into findcut2",13,k,1)
rmax = rmaxi
rmin = rmini
30 rtest = (rmax + rmin)/2
call checklev(f,g,nbin,mc,ml,mr,k,idir,rtest,a1,a2)
if (diff .lt. 1D-50) then
write (*,*) 'diff = ',diff,'. Div by 0 in findcut2.'
end if
if ((abs(a2-diff)/diff) .lt. 0.0001) goto 500
if (diff .gt. a2) then
rmin = rtest
else
rmax = rtest
end if
goto 30
500 return
end
c April 16, 1992
c find level to cut (one side) equalize mass above and below cut
c
subroutine findcut(f,g,nbin,mc,ml,mr,k,idir,rmaxi,rmini,rtest)
double precision f(nbin), g(nbin), rmaxi, rmini, rtest, rmax,
+ rmin, a1, a2
integer mc(100), ml(100), mr(100), nbin, k, idir
c call intpr("Into findcut",12,k,1)
rmax = rmaxi
rmin = rmini
30 rtest = (rmax + rmin)/2
call checklev(f,g,nbin,mc,ml,mr,k,idir,rtest,a1,a2)
if (a1 .lt. 1D-50) then
write (*,*) 'a1 = ',a1,'. Div by 0 in findcut.'
end if
if ((abs(a1-a2)/a1) .lt. 0.0001) goto 500
if (a1 .gt. a2) then
rmin = rtest
else
rmax = rtest
end if
goto 30
500 return
end
c April 16, 1992
c get values above and below level
subroutine checklev(f,g,nbin,mc,ml,mr,k,idir,rtest,a1,a2)
double precision f(nbin), g(nbin), rtest, a1, a2, aid
integer mc(100), ml(100), mr(100), nbin, k, idir, ml1, mr1, ml2,
+ mr2, i
c call dblepr("Into checklev",13,rtest,1)
call level(f, rtest, nbin, ml(k), mr(k), mc(k), ml1, mr1)
call level(f, rtest, nbin, mc(k+min(0,idir)),mc(k+max(0,idir)),
+ mr(k+min(0,idir)), ml2,mr2)
c call intpr("ml1",3,ml1,1)
c call intpr("mr1",3,mr1,1)
c call intpr("ml2",3,ml2,1)
c call intpr("mr2",3,mr2,1)
do 40 i = 1, nbin
g(i) = f(i)
40 continue
do 50 i = min(ml1,ml2), max(mr1, mr2)
g(i) = rtest
50 continue
a1 = aid(f,g,ml1,mr1,nbin)
a2 = aid(f,g,ml2,mr2,nbin)
return
end
subroutine level(f, rlev, nbin, mlin, mrin, mcin, mlout, mrout)
double precision f(nbin), rlev
integer nbin, mlin, mrin, mcin, mlout, mrout, i
c call intpr("Into level",10,mlin,1)
if (f(mcin) .gt. f(mlin)) then
i = mlin + 1
10 if (f(i) .ge. rlev) goto 20
i = i + 1
goto 10
20 continue
mlout = i
30 i = i + 1
if (f(i) .le. rlev) goto 40
goto 30
40 continue
mrout = i - 1
else
i = mlin + 1
50 if (f(i) .le. rlev) goto 60
i = i + 1
goto 50
60 continue
mlout = i
70 i = i + 1
if (f(i) .ge. rlev) goto 80
goto 70
80 continue
mrout = i - 1
endif
return
end