## A knapsack riddle [#2]?

Posted in Kids, R, Statistics with tags , , , on February 17, 2017 by xi'an

Still about this allocation riddle of the past week, and still with my confusion about the phrasing of the puzzle, when looking at a probabilistic interpretation of the game, rather than for a given adversary’s y, the problem turns out to search for the maximum of

$\mathbb{E}[L(x,Y)]=\sum_{i=1}^{10} i\{P(Y_ix_i)\}$

where the Y’s are Binomial B(100,p). Given those p’s, this function of x is available in closed form and can thus maximised by a simulated annealing procedure, coded as

utility=function(x,p){
ute=2*pbinom(x[1]-1,100,prob=p[1])+
dbinom(x[1],100,p[1])
for (i in 2:10)
ute=ute+2*i*pbinom(x[i]-1,100,prob=p[i])+
i*dbinom(x[i],100,p[i])
return(ute)}
#basic term in utility
baz=function(i,x,p){
return(i*dbinom(x[i],100,p[i])+
i*dbinom(x[i]-1,100,p[i]))}
#relies on a given or estimated p
x=rmultinom(n=1,siz=100,prob=p)
maxloz=loss=0
newloss=losref=utility(x,p)
#random search
T=1e3
Te=1e2
baza=rep(0,10)
t=1
while ((t<T)||(newloss>loss)){
loss=newloss
i=sample(1:10,1,prob=(10:1)*(x>0))
#moving all other counters by one
xp=x+1;xp[i]=x[i]
#corresponding utility change
for (j in 1:10) baza[j]=baz(j,xp,p)
proz=exp(log(t)*(baza-baza[i])/Te)
#soft annealing move
j=sample(1:10,1,prob=proz)
if (i!=j){ x[i]=x[i]-1;x[j]=x[j]+1}
newloss=loss+baza[j]-baza[i]
if (newloss>maxloz){
maxloz=newloss;argz=x}
t=t+1
if ((t>T-10)&(newloss<losref)){
t=1;loss=0
x=rmultinom(n=1,siz=100,prob=p)
newloss=losref=utility(x,p)}}


which seems to work, albeit not always returning the same utility. For instance,

> p=cy/sum(cy)
> utility(argz,p)
[1] 78.02
> utility(cy,p)
[1] 57.89


or

> p=sy/sum(sy)
> utility(argz,p)
[1] 82.04
> utility(sy,p)
[1] 57.78


Of course, this does not answer the question as intended and reworking the code to that purpose is not worth the time!

## Le Monde puzzle [#958]

Posted in Books, Kids, R with tags , , , on April 11, 2016 by xi'an

A knapsack Le Monde mathematical puzzle:

Given n packages weighting each at most 5.8kg for a total weight of 300kg, is it always possible to allocate these packages  to 12 separate boxes weighting at most 30kg each? weighting at most 29kg each?

This can be checked by brute force using the following R code

#generate packages
paca=runif(1,0,5.8)
while (sum(paca)<300){
paca=c(paca,runif(1,0,5.8))}
paca=paca[-length(paca)]
paca=c(paca,300-sum(paca))


and

#check if they can be partitioned into
#12 groups of weight less than 30
box=vector(mode="list",length=12)
#random allocation
alloc=sample(1:12,length(paca),rep=TRUE)
content=rep(0,12)
for (i in 1:12){
box[[i]]=paca[alloc==i]
content[i]=sum(box[[i]])}
content=content*(content>0)
#wrong allocation
while (max(content)>30){
i=sample(1:12,1,prob=content)
j=sample(1:length(box[[i]]),1,prob=box[[i]])
#reallocation
k=sample(1:12,1,prob=max(content)-content)
while (k==i){
k=sample(1:12,1,prob=max(content)-content)}
content[i]=content[i]-box[[i]][j]
content[i]=content[i]*(content[i]>0)
content[k]=content[k]+box[[i]][j]
box[[k]]=c(box[[k]],box[[i]][j])
box[[i]]=box[[i]][-j]}


repeatedly and could not find an endless while loop. (Empty boxes sometimes lead to negative contents, hence my fix setting negative contents to zero.) But neither did I find an issue when the upper bound on the weight is 29kg… So it is either possible or rarely impossible! The R code immediately gets stuck when setting the bound at 25kg.

After reading the solution of April 13 in Le Monde, it appears that there is a counter example for the 29kg limit, namely 60 packages weighting 4.91kg plus one package weighting 5.4kg, since the first 60 packages fit inside 12 boxes and the last one is left out.