# 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 ri and and compares the minimum produced by the code below to the global search over all .for the minimum of .

#  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    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)
}
}
break 1
low = mht+1
large = mlt-1
if (low>large)
break 1
gn = gnt+gmix
gp = gpt+gplx
ml = mlt
}
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