c Some sample Fortran code for the binmedian algorithm. c Note: I haven't written the code for n even yet. real function bmedian(n,x) integer n real x(n) integer i,bottom,counts(0:1000),bin,k,count,medbin,r,oldbin,j real mu,sigma,scalefac,leftend,rghtend,a,oldscale,oldleft logical samepts; c Compute the mean and standard deviation mu = sum(x)/n sigma = sqrt(dot_product(x-mu,x-mu)/n) c Bin x across the interval [mu-sigma, mu+sigma] bottom = 0 counts = 0 scalefac = 1000/(2*sigma) leftend = mu-sigma rghtend = mu+sigma do 1 i = 1,n if (x(i).lt.leftend) then bottom = bottom+1 else if (x(i).lt.rghtend) then bin = int((x(i)-leftend) * scalefac) counts(bin) = counts(bin)+1 endif 1 continue c If n is odd if (mod(n,2).eq.1) then k = (n+1)/2 r = 1 c Beginning of recursive step: c Find the bin that contains the median, and the order c of the median within that bin 2 count = bottom do 3 i = 0,1000 count = count+counts(i) if (count.ge.k) then medbin = i k = k - (count-counts(i)) goto 4 endif 3 continue 4 bottom = 0 counts = 0 oldscale = scalefac oldleft = leftend scalefac = 1000*oldscale leftend = medbin/oldscale + oldleft rghtend = (medbin+1)/oldscale + oldleft c Determine which points map to medbin, and put c them in spots r,...n i = r r = n+1 5 if (i.lt.r) then oldbin = int((x(i)-oldleft) * oldscale) if (oldbin.eq.medbin) then r = r-1 c Swap x(i) and x(r) a = x(i) x(i) = x(r) x(r) = a c Re-bin on a finer scale if (x(r).lt.leftend) then bottom = bottom+1 else if (x(r).lt.rghtend) then bin = int((x(r)-leftend) * scalefac) counts(bin) = counts(bin)+1 endif else i = i+1 endif goto 5 endif c Stop if all points in medbin are the same samepts = .TRUE. do 6 i = r+1,n if (x(i).ne.x(r)) then samepts = .FALSE. goto 7 endif 6 continue 7 if (samepts) then bmedian = x(r) return endif c Stop if there's <= 20 points left if ((n+1-r).le.20) then goto 8 endif goto 2 c Perform insertion sort on the remaining points, c and then pick the kth smallest 8 do 9 i = r+1,n a = x(i) do 10 j = i-1,r,-1 if (x(j).lt.a) goto 11 x(j+1) = x(j) 10 continue j = r-1 11 x(j+1) = a 9 continue bmedian = x(r-1+k) return c If n is even else bmedian = 0 return endif end