CODE FOR WEIGHTED MEDIAN

I adapted the code below from FGDP to find the weighted median. Theoretically this code requires about three times as many arithmetical operators as a weighted mean. Instead of N multiplications and additions, the weighted median is expected to require about 3N multiplications, subtractions, and polarity tests, as well as some lessor number of element interchanges.

If you have the electronic document you can test the code by typing gmake minabs.test which uses random numbers for components r_i and $\Delta r_i$and compares the minimum produced by the code below to the global search over all $\alpha_j=r_j/\Delta r_j$.for the minimum of $\sum_i\vert r_i+\alpha_j\Delta r_i\vert$.

#  minimize(alpha)  =   sum  abs( f(i) + w(i)*alpha )
#                        i
subroutine minabs( n, f, w, alpha)
integer i,n,low,large,ml,mh,mlt,mht,itry,j
real gn,gp,t,alpha,er,gnt,gpt,gplx,gmix,small,grad
          real    w(n),f(n)
temporary real    g(n)
temporary integer k(n)
small = 0.
do i=1,n
	k(i) = i
low = 1
large = n
ml = n
mh = 1
gn = 0.
gp = 0.
do itry = 1,n {
	j = k( low+mod( (large-low)/3+itry, large-low+1) )
	if( abs(w(j)) != 0.) {
		t = f(j)/(w(j))
		f(j) = w(j)*t
		do i = low,large {
			j = k(i)
			er = f(j)-w(j)*t
			g(j) = 0.
			if     ( er>small)
					g(j) = -w(j)
			else if( er<(-small))
					g(j) = +w(j)
			else
					g(j) = 0.
			}
		call split(low,large,k,g,mlt,mht)
		gnt = gn
		do i = low,mlt
			gnt = gnt+g(k(i))
		gpt = gp
		do i = mht,large
			gpt = gpt+g(k(i))
		gplx = 0.
		gmix = 0.
		do i = mlt,mht {
			j = k(i)
			if (w(j)<0.) {
				gplx = gplx-w(j)
				gmix = gmix+w(j)
				}
			if (w(j)>0.) {
				gplx = gplx+w(j)
				gmix = gmix-w(j)
				}
			}
		grad = gnt+gpt
		if ((grad+gplx)*(grad+gmix)<0.)
				break 1
		if (grad>=0.)
				low = mht+1
		if (grad<=0)
				large = mlt-1
		if (low>large)
				break 1
		if (grad>=0.)
				gn = gnt+gmix
		if (grad<=0.)
				gp = gpt+gplx
		if (grad+gplx==0.)
				ml = mlt
		}
	if (grad+gmix==0.)
			mh = mht
	}
alpha = - t
return; end

subroutine split(low,large,k,g,ml,mh)
#	given g(k(i)),i=low,large
#	then rearrange k(i),i=low,large and find ml and mh so that
#	(g(k(i)),i=low,(ml-1)) < 0. and
#	(g(k(i)),i=ml,mh)=0. and
#	(g(k(i)),i=(mh+1),large) > 0.
real    g(1)
integer k(1)
integer low,large,ml,mh,keep,i,ii
ml = low
mh = large
repeat {
	ml = ml-1
	repeat
		ml = ml+1
		until(g(k(ml))>=0)
	repeat {
		mh = mh+1
		repeat
			mh = mh-1
			until(g(k(mh))<=0)
		keep = k(mh)
		k(mh) = k(ml)
		k(ml) = keep
		if (g(k(ml))!=g(k(mh)))
			break 1		# Break out of enclosing "repeat{"
		do i = ml,mh {
			ii = i
			if (g(k(i))!=0.0)
				go to 30
			}
		break 2		# Break out of two enclosing "repeat{"'s
		30  keep = k(mh)
		k(mh) = k(ii)
		k(ii) = keep
		}
	}
return
end