From sujith@math.rutgers.edu  Sat Dec 18 01:10:40 2004
Received: from localhost (sujith@localhost)
	by math.rutgers.edu (8.11.7p1+Sun/8.8.8) with ESMTP id iBI6Ad525482
	for <zeilberg@math.rutgers.edu>; Sat, 18 Dec 2004 01:10:39 -0500 (EST)
Date: Sat, 18 Dec 2004 01:10:39 -0500 (EST)
From: Sujith Vijay <sujith@math.rutgers.edu>
To: Doron Zeilberger <zeilberg@math.rutgers.edu>
Subject: Knock 'em Down
Message-ID: <Pine.GSO.4.44.0412180039190.22836-100000@math.rutgers.edu>
MIME-Version: 1.0
Content-Type: TEXT/PLAIN; charset=US-ASCII
Status: RO
Content-Length: 4094

Dear Professor Zeilberger,

Here is what I have:

I have a heuristic, tested only for three piles, but can be generalised,
which computes near-optimal allocations quickly. Brute force takes way too
long, the idea is to reduce it to a problem of two piles.

So suppose we have three bins, with probabilities p = 1/7, q = 2/7 and
1-p-q = 4/7 and we have to allocate 50 tokens.

We could search over all possible allocations, but this is time consuming.

Instead, we consider the first pile against the other two together, and
pretend that we are in a two pile situation with probabilities 1/7 and
6/7. This calls for 5 tokens in the first pile.

Now we consider the second pile against the other two, with probabilities
2/7 and 5/7. This would mean 12 tokens in the second pile.

Finally, we consider the third pile against the first two, and put 29
tokens in the third pile.

That still leaves 4 tokens, but we appeal to recursion. It turns out that
(0,1,3) is the best allocation for 4 tokens, given the above
probabilities.  Adding this to the pile, we get the allocation (5,13,32).

Miraculously, we never overshoot (actually, it isn't that big a miracle,
at most one of the probabilities can be more than 1/2, and optimal
allocations for two bins can be shown to drift in opposite directions from
the `mean' allocation, so most bins go below the mean).

The optimal allocation is (5,14,31), but this took 29 seconds of
computation time, while the heuristic results were instantaneous. I have
tried more than thirty examples, with various probabilities, and up to 70
tokens, and the R.M.S. error over the three co-ordinates has always been
sqrt(2) or less, the worst case error being (1,1,-2).

I have to appear for the foreign language exam next week, but if I clear
that, I hope to take another look at these results. I have already spent a
few days trying to prove error bounds, though nothing has been
forthcoming.

Here is the relevant Maple Code:

-----

Emp:=proc(m,p):
if m=0 then
  0:
else
 (p*Emp(m-1,p)+1)/p:
fi:
end:


Emnp:=proc(m::nonnegint,n::nonnegint,p)
option remember:

if m>0 and n>0 then
 p*Emnp(m-1,n,p) +(1-p)* Emnp(m,n-1,p)+1:
elif m>0 and n=0 then
Emp(m,p):
elif m=0 and n>0 then
Emp(n,1-p):
else
0:
fi:
end:


E3:=proc(x::nonnegint,y::nonnegint,z::nonnegint,p,q)
option remember:
if x>0 and y>0 and z > 0 then
p*E3(x-1,y,z,p,q)+q*E3(x,y-1,z,p,q)+(1-p-q)*E3(x,y,z-1,p,q)+1:
elif x=0 and y>0 and z>0 then
(q*E3(0,y-1,z,p,q)+(1-p-q)*E3(0,y,z-1,p,q)+1)/(1-p):
elif y=0 and x>0 and z>0 then
(p*E3(x-1,0,z,p,q)+(1-p-q)*E3(x,0,z-1,p,q)+1)/(1-q):
elif z=0 and x>0 and y>0 then
(p*E3(x-1,y,0,p,q)+q*E3(x,y-1,0,p,q)+1)/(p+q):
elif x=0 and y=0 and z>0 then z/(1-p-q):
elif x=0 and z=0 and y>0 then y/q:
elif y=0 and z=0 and x>0 then x/p:
else 0:
fi:
end:

Opt3:=proc(n::nonnegint,p,q) local i,j,k,besti,bestj,best,val:
besti:=0:
bestj:=0:
best:=E3(0,0,n,p,q):
for i from 0 to n do
for j from 0 to n do
k:=n-i-j:
if k >= 0 then val:=E3(i,j,k,p,q): fi:
   if val < best then
      besti:=i:
      bestj:=j:
      best:=val:
   fi:
od:
od:
[besti,bestj,evalf(best)]:
end:

listsub:=proc(ls1,ls2) local i:
[seq(ls1[i]-ls2[i],i=1..min(nops(ls1),nops(ls2)))]:
end:

listadd:=proc(ls1,ls2) local i:
[seq(ls1[i]+ls2[i],i=1..min(nops(ls1),nops(ls2)))]:
end:

Opt2:=proc(n::nonnegint,p) local i,j,besti,best,val:
besti:=0:
best:=Emnp(0,n,p):
for i from 0 to n do
j:=n-i:
if j >= 0 then val:=Emnp(i,j,p): fi:
   if val < best then
      besti:=i:
      best:=val:
   fi:
od:
[besti,evalf(best)]:
end:

Heur:=proc(n::nonnegint,p,q) local a,b,c,count,k,s,ls,ls1,ls2,res,sum:
k:=n:
ls2:=[0,0,0]:
res:=k:
while res > 2 do
a:=Opt2(res,p)[1]:
b:=Opt2(res,q)[1]:
c:=Opt2(res,1-p-q)[1]:
ls2:=listadd(ls2,[a,b,c]):
sum:=a+b+c:
res:=res-sum:
if res < 0 then return FAIL:fi:
od:
ls2:=listadd(ls2,Opt3(res,p,q)):
ls2[3]:=evalf(E3(ls2[1],ls2[2],n-ls2[1]-ls2[2],p,q)):
ls2:
end:

-----