Archive for R

Le Monde puzzle [#1021]

Posted in Books, Kids, R with tags , , , , , on September 17, 2017 by xi'an

A puzzling Le Monde mathematical puzzle for which I could find no answer in the allotted time!:

A most democratic electoral system allows every voter to have at least one representative by having each of the N voters picking exactly m candidates among the M running candidates and setting the size n of the representative council towards this goal, prior to the votes. If there are M=25 candidates, m=10 choices made by the voters, and n=10 representatives, what is the maximal possible value of N? And if N=55,555 and M=33, what is the minimum value of n for which m=n is always possible?

I tried a brute force approach by simulating votes from N voters at random and attempting to find the minimal number of councillors for this vote, which only provides an upper bound of the minimum [for one vote], and a lower bound in the end [over all votes]. Something like

for (i in 1:N) votz[i,]=sample(1:M,n)
#exploration by majority
  remz=1:N;conz=NULL
  while (length(remz)>0){
    seatz=order(-hist(votz[remz,],
    breaks=(0:M)+0.5,plot=FALSE)$density)[1]
    conz=c(conz,seatz);nuremz=NULL
    for (v in remz)
      if (!(seatz%in%votz[v,])) nuremz=c(nuremz,v)
    remz=nuremz}
  solz=length(conz)
#exploration at random
   kandz=matrix(0,N,M)
   for (i in 1:N) kandz[i,votz[i,]]=1
   for (t in 1:1e3){
#random choice of councillors
    zz=sample(c(0,1),M,rep=TRUE)
    while (min(kandz%*%zz)!=1)
      zz=sample(c(0,1),M,rep=TRUE)
    solz=min(solz,sum(zz))
#random choice of remaining councillor per voter
    remz=1:N;conz=NULL
    while (length(remz)>0){
      seatz=sample(votz[remz[1],],1)
      conz=c(conz,seatz);nuremz=NULL
      for (i in remz)
        if (!(seatz%in%votz[i,])) nuremz=c(nuremz,i)
      remz=nuremz}
    solz=min(solz,length(conz))}
maxz=max(solz,maxz)}

which leads to a value near N=4050 for the first question, with 0% confidence… Obviously, the problem can be rephrased as a binary integer linear programming problem of the form

n= \max_A \min_{c;\,\min Ac=1}\mathbf{1}^\text{T}c

where A is the NxM matrix of votes and c is the vector of selected councillors. But I do not see a quick way to fix it!

Le Monde puzzle [#1019]

Posted in Books, Kids with tags , , , , , , on September 7, 2017 by xi'an

A gamey (and verbose) Le Monde mathematical puzzle:

A two-player game involves n+2 cards in a row, blue on one side and red on the other. Each player can pick any blue card among the n first ones and flip it plus both following ones. The game stops when no blue card is left to turn. The gain for the last player turning cards is 20-t, where t is the number of times cards were flipped, with gain t for its opponent. Both players aim at maximising their gain.

1. When n=4 and all cards are blue, can the first player win? If not, what is the best score for this player?

2. Among all 16 configurations at start, how many lead to the first player to win?

3. When n=10 and all cards are blue, how many cards are flipped an odd number of times for the winning configuration?

The first two questions can easily be processed by an R code like the following recursive function:

liplop <- function(x,n,i){
  if (max(x[1:n])==0){
    return(i)
  }else{
    sol=NULL
    for (j in (1:n)[x[1:n]==1]){
      y=x;y[j:(j+2)]=1-y[j:(j+2)]
      sol=c(sol,20-liplop(y,n,i+1))}
    return(max(sol))}}

Returning

> liplop(rep(1,6),4,0)
[1] 6

Meaning the first player cannot win, by running at most six rounds. Calling the same function for all 4⁴=16 possible configurations leads to 8 winning ones:

[1] 0 0 0 1
[1] 0 0 1 1
[1] 0 1 0 1
[1] 0 1 1 1
[1] 1 0 0 0
[1] 1 0 1 0
[1] 1 1 0 0
[1] 1 1 1 0

Solving the same problem with n=10 is not feasible with this function. (Even n=6 seems out of reach!)

Le Monde puzzle [#1018]

Posted in Books, Kids, R with tags , , , , , on August 29, 2017 by xi'an

An arithmetic Le Monde mathematical puzzle (that first did not seem to involve R programming because of the large number of digits in the quantity involved):

An integer x with less than 100 digits is such that adding the digit 1 on both sides of x produces the integer 99x.  What are the last nine digits of x? And what are the possible numbers of digits of x?

The integer x satisfies the identity

10^{\omega+2}+10x+1=99x

where ω is the number of digits of x. This amounts to

10….01 = 89 x,

where there are ω zeros. Working with long integers in R could bring an immediate solution, but I went for a pedestrian version, handling each digit at a time and starting from the final one which is necessarily 9:

#multiply by 9
rap=0;row=NULL
for (i in length(x):1){
prud=rap+x[i]*9
row=c(prud%%10,row)
rap=prud%/%10}
row=c(rap,row)
#multiply by 80
rep=raw=0
for (i in length(x):1){
prud=rep+x[i]*8
raw=c(prud%%10,raw)
rep=prud%/%10}
#find next digit
y=(row[1]+raw[1]+(length(x)>1))%%10

returning

7 9 7 7 5 2 8 0 9

as the (only) last digits of x. The same code can be exploited to check that the complete multiplication produces a number of the form 10….01, hence to deduce that the length of x is either 21 or 65, with solutions

[1] 1 1 2 3 5 9 5 5 0 5 6 1 7 9 7 7 5 2 8 0 9
[1] 1 1 2 3 5 9 5 5 0 5 6 1 7 9 7 7 5 2 8 0 8 9 8 8 7 6 4 0 4 4 9 4 3 8 2 0 2 2
[39] 4 7 1 9 1 0 1 1 2 3 5 9 5 5 0 5 6 1 7 9 7 7 5 2 8 0 9

The maths question behind is to figure out the powers k of 10 such that

10^k\equiv -1 \text{ mod } (89)

For instance, 10²≡11 mod (89) and 11¹¹≡88 mod (89) leads to the first solution ω=21. And then, since 10⁴⁴≡1 mod (89), ω=21+44=65 is another solution…

Le Monde puzzle [#1707]

Posted in Books, Kids, R with tags , , , , , on July 28, 2017 by xi'an

A geometric Le Monde mathematical puzzle:

  1. Given a pizza of diameter 20cm, what is the way to cut it by two perpendicular lines through a point distant 5cm from the centre towards maximising the surface of two opposite slices?
  2.  Using the same point as the tip of the four slices, what is the way to make four slices with equal arcs in four cuts from the tip again towards maximising the surface of two opposite slices?

For both questions, I did not bother with the maths but went itself to a discretisation of the disk, counting the proportion of points within two opposite slices and letting the inclination of these slices move from zero to π/2. Unsurprisingly, for the first question, the answer is π/4, given that there is no difference between both surfaces at angles 0 and π/2. My R code is as follows, using (5,0) as the tip:

M=100
surfaz=function(alpha){
surfz=0
cosal=cos(alpha);sinal=sin(alpha)
X=Y=seq(-10,10,le=M)
Xcosal=(X-5)*cosal
Xsinal=(X-5)*sinal
for (i in 1:M){
norm=sqrt(X[i]^2+Y^2)
scal1=Xsinal[i]+Y*cosal
scal2=-Xcosal[i]+Y*sinal
surfz=surfz+sum((norm<=10)*(scal1*scal2>0))}
return(4*surfz/M/M/pi)}

The second puzzle can be solved by a similar code, except that the slice area between two lines has to be determined by a cross product:

surfoz=function(alpha,ploz=FALSE){
  sinal=sin(alpha);cosal=cos(alpha)
  X=Y=seq(-10,10,le=M)
  frsterm=cosal*(10*cosal-5)+sinal*(10*sinal-5)
  trdterm=cosal*(10*cosal+5)+sinal*(10*sinal+5)
  surfz=0
  for (i in 1:M){
    norm=sqrt(X[i]^2+Y^2)
    scal1=(10*(Y[i]-5)*cosal-(10*sinal-5)*X)*frsterm
    scal2=-(-10*(Y[i]-5)*sinal-(10*cosal-5)*X)*frsterm
    scal3=(-10*(Y[i]-5)*cosal+(10*sinal+5)*X)*trdterm
    scal4=-(10*(Y[i]-5)*sinal+(10*cosal+5)*X)*trdterm
    surfz=surfz+sum((norm<=10)* 
    ((scal1>0)*(scal2>0)+
     (scal3>0)*(scal4>0)))}
 return(4*surfz/M/M/pi)}

a code that shows that all cuts lead to identical surfaces for bot sets of slices. A fairly surprising result!

 

Le Monde puzzle [#1016]

Posted in Books, Kids with tags , , , on July 16, 2017 by xi'an

An even more straightforward Le Monde mathematical puzzle that took a few minutes to code in the train to Cambridge:

  1. Breaking {1,…,8} into two sets of four integrals, what is (or are) the division into two groups of equal size such that the sums of the squared terms from each are equal? Same question for the set {21,…,28}.
  2.  Considering the integers from 1 to 12, how many divisions into two groups of size six satisfy the above property? Same question when the two groups are of different sizes.

The first code is

nop=TRUE
while (nop){
 s=sample(1:8)
 nop=(sum(s[1:4]^2)!=sum(s[5:8]^2))}

with result

1 6 4 7

while the second set leads to the unique [drifted] solution (up to symmetries)

21 24 26 27

and the divisions for the larger set {1,…,12} is unique in the equal case, and are four in the unequal case.

RNG impact on MCMC [or lack thereof]

Posted in Books, R, Statistics, Travel, University life with tags , , , , , , , on July 13, 2017 by xi'an

Following the talk at MCM 2017 about the strange impact of the random generator on the outcome of an MCMC generator, I tried in Montréal airport the following code on the banana target of Haario et al. (1999), copied from Soetaert and Laine and using the MCMC function of the FME package:

library(FME)
Banana <- function (x1, x2) {
 return(x2 - (x1^2+1)) }
pmultinorm <- function(vec, mean, Cov) {
 diff <- vec - mean
 ex <- -0.5*t(diff) %*% solve(Cov) %*% diff
 rdet <- sqrt(det(Cov))
 power <- -length(diff)*0.5
 return((2.*pi)^power / rdet * exp(ex)) }
BananaSS <- function (p) {
 P <- c(p[1], Banana(p[1], p[2]))
 Cov <- matrix(nr = 2, data = c(1, 0.9, 0.9, 1))
N=1e3
ejd=matrix(0,4,N)
RNGkind("Mars")
for (t in 1:N){
  MCMC <- modMCMC(f = BananaSS, p = c(0, 0.7), 
  jump = diag(nrow = 2, x = 5), niter = 1e3)
  ejd[1,t]=mean((MCMC$pars[-1,2]-MCMC$pars[1,2])^2)}

since this divergence from the initial condition seemed to reflect the experiment of the speaker at MCM 2017. Unsurprisingly, no difference came from using the different RNGs in R (which may fail to contain those incriminated by the study)…

easy riddle

Posted in Books, Kids, R with tags , , , , , on July 12, 2017 by xi'an

From the current Riddler, a problem that only requires a few lines of code and a few seconds of reasoning. Or not.

N households each stole the earnings from one of the (N-1) other households, one at a time. What is the probability that a given household is not burglarised? And what are the expected final earnings of each household in the list, assuming they all start with $1?

The first question is close to Feller’s enveloppe problem in that

\left(1-\frac{1}{N-1}\right)^{N-1}

is close to exp(-1) for N large. The second question can easily be solved by an R code like

N=1e3;M=1e6
fina=rep(1,N)
for (v in 1:M){
 ordre=sample(1:N)
 vole=sample(1:N,N,rep=TRUE)
 while (min(abs(vole-(1:N)))==0)
  vole[abs(vole-(1:N))==0]=sample(1:N,
     sum(vole-(1:N)==0))
 cash=rep(1,N)
 for (t in 1:N){
  cash[ordre[t]]=cash[ordre[t]]+cash[vole[t]];cash[vole[t]]=0}
 fina=fina+cash[ordre]}

which returns a pretty regular exponential-like curve, although I cannot figure the exact curve beyond the third burglary. The published solution gives the curve

{\frac{N-2}{N-1}}^{999}\times 2+{\frac{1}{N-1}}^{t-1}\times{\frac{N-1}{N}}^{N-t}\times\frac{N}{N-1}

corresponding to the probability of never being robbed (and getting on average an extra unit from the robbery) and of being robbed only before robbing someone else (with average wealth N/(N-1)).