Archive for optimisation

oxwasp@amazon.de

Posted in Books, Kids, pictures, Running, Statistics, Travel, University life with tags , , , , , , , , , , , , , , , , , on April 12, 2017 by xi'an

The reason for my short visit to Berlin last week was an OxWaSP (Oxford and Warwick Statistics Program) workshop hosted by Amazon Berlin with talks between statistics and machine learning, plus posters from our second year students. While the workshop was quite intense, I enjoyed very much the atmosphere and the variety of talks there. (Just sorry that I left too early to enjoy the social programme at a local brewery, Brauhaus Lemke, and the natural history museum. But still managed nice runs east and west!) One thing I found most interesting (if obvious in retrospect) was the different focus of academic and production talks, where the later do not aim at a full generality or at a guaranteed improvement over the existing, provided the new methodology provides a gain in efficiency over the existing.

This connected nicely with me reading several Nature articles on quantum computing during that trip,  where researchers from Google predict commercial products appearing in the coming five years, even though the technology is far from perfect and the outcome qubit error prone. Among the examples they provided, quantum simulation (not meaning what I consider to be simulation!), quantum optimisation (as a way to overcome multimodality), and quantum sampling (targeting given probability distributions). I find the inclusion of the latest puzzling in that simulation (in that sense) shows very little tolerance for errors, especially systematic bias. It may be that specific quantum architectures can be designed for specific probability distributions, just like some are already conceived for optimisation. (It may even be the case that quantum solutions are (just next to) available for intractable constants as in Ising or Potts models!)

Le Monde puzzle [#1002]

Posted in Kids, R with tags , , , , , , , on April 4, 2017 by xi'an

For once and only because it is part of this competition, a geometric Le Monde mathematical puzzle:

Given both diagonals of lengths p=105 and q=116, what is the parallelogram with the largest area? and when the perimeter is furthermore constrained to be L=290?

This made me jump right away to the quadrilateral page on Wikipedia, which reminds us that the largest area occurs when the diagonals are orthogonal, in which case it is A=½pq. Only the angle between the diagonals matters. Imposing the perimeter 2s in addition is not solved there, so I wrote an R code looking at all the integer solutions, based on one of the numerous formulae for the area, like ½pq sin(θ), where θ is the angle between both diagonals, and discretising in terms of the fractions of both diagonals at the intersection, and of the angle θ:

p=105
q=116
s=145
for (alpha in (1:500)/1000){
 ap=alpha*p;ap2=ap^2;omap=p-ap;omap2=omap^2
 for (beta in (1:999)/1000){
  bq=beta*q;bq2=bq^2;ombq=q-bq;ombq2=ombq^2
  for (teta in (1:9999)*pi/10000){
   d=sqrt(ap2+bq2-2*ap*bq*cos(teta))
   a=sqrt(ap2+ombq2+2*ap*ombq*cos(teta))
   b=sqrt(omap2+ombq2-2*omap*ombq*cos(teta))
   c=sqrt(omap2+bq2+2*omap*bq*cos(teta))
   if (abs(a+b+c+d-2*s)<.01){
     if (p*q*sin(teta)<2*maxur){
       maxur=p*q*sin(teta)/2
       sole=c(a,b,c,d,alpha,beta,teta)}}}}

This code returned an area of 4350, to compare with the optimal 6090 (which is recovered by the above R code when the diagonal lengths are identical and the perimeter is the one of the associated square). (As Jean-Louis Foulley pointed out to me, this area can be found directly by assuming the quadrilateral is a parallelogram and maximising in the length of one side.)

A knapsack riddle [#2]?

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

gear

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_i<x_i)-P(Y_i>x_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!

a knapsack riddle?

Posted in Books, pictures, R, Statistics, Travel with tags , , , , , , on February 13, 2017 by xi'an

gear

The [then current now past] riddle of the week is a sort of multiarmed bandits optimisation. Of sorts. Or rather a generalised knapsack problem. The question is about optimising the allocation of 100 undistinguishable units to 10 distinct boxes against a similarly endowed adversary, when the loss function is

L(x,y)=(x_1>y_1)-(x_1<y_1)+...+10((x_{10}>y_{10})-(x_{10}<y_{10}))

and the distribution q of the adversary is unknown. As usual (!), the phrasing of the riddle is somewhat ambiguous but I am under the impression that the game is played sequentially, hence that one can learn about the distribution of the adversary, at least when assuming this adversary keeps the same distribution q at all times. Continue reading

optimality stands in the eye of the riddler

Posted in Books, Kids, pictures, Statistics with tags , , , , , , on March 22, 2016 by xi'an

When looking at US elections on FiveThirtyEight, I came across this riddle:

Two players go (…) into separate booths (…) and a random number between zero and one appears on a screen (…) chosen from a standard uniform distribution. They can choose to keep that first number, or to (…) get a second random number, which they must keep (…) The prize (…) is awarded to the player who kept the higher number (…) Which number is the optimal cutoff for players to discard their first number and choose another? [The Riddler, Mar 4, 2016]

While the solution is now available, I wonder at the wording of this riddle, where “optimality” is not spelled out. Unless I missed some key entry, as it often happens with puzzles… Assuming both players use the same “optimal” cut-off C to make their decision to keep or drop the original uniform, the probability that one does better than the other is exactly ½ since both outcomes are iid. When considering the expected value of the number kept by a player, a simple integral shows that this value is

½(1-C²+C),

which is maximal for C=½. If one considers instead the median of the number Z kept by a player, a bit more computation leads to the function med(Z) = 1/2C when C²>1/2 and (½+C)/(1+C) when C²<1/2, which is maximal for C²=½.

“…using the golden ratio gives the best chance to win the gold bullion!” [The Riddler, Mar 4, 2016]

Incidentally (or not), the solution on 538 is also difficult to understand in that it starts with the event that the first draw is C, which is an event of probability zero. However, when trying to optimise the choice of C for one player, given that the other player has a known cuttoff of D, I found that the only case when C=D coincides with the solution proposed on Riddler, namely ½(√5-1). To check whether or not this derivation was correct, I also plotted the theoretical (right) versus the empirical (left) versions of the winning probability:

riddlerThere is no difference between the two. But the diagonal is exactly at the .5 level, as expected:riddlerinthestormwith an interesting second curve at the .5 probability level. These two level sets happen to meet at the “golden number” solution, ½(√5-1), which comes as a confirmation of my earlier remark. [Another question connected with this riddle would be to figure out the value of D used by the other player from a sequence of games.]

Approximate Maximum Likelihood Estimation

Posted in Books, Mountains, pictures, Statistics, Travel, University life with tags , , , , , , , , , on September 21, 2015 by xi'an

linz3Bertl et al. arXived last July a paper on a maximum likelihood estimator based on an alternative to ABC techniques. And to indirect inference. (One of the authors in et al. is Andreas Futschik whom I visited last year in Linz.) Paper that I only spotted when gathering references for a reading list on ABC… The method is related to the “original ABC paper” of Diggle and Gratton (1984) which, parallel to Rubin (1984), contains in retrospect the idea of ABC methods. The starting point is stochastic approximation, namely the optimisation of a function of a parameter θ when written as an expectation of a random variable Y, E[Y|θ], as in the Kiefer-Wolfowitz algorithm. However, in the case of the likelihood function, there is rarely an unbiased estimator and the authors propose instead to use a kernel density estimator of the density of the summary statistic. This means that, at each iteration of the Kiefer-Wolfowitz algorithm, two sets of observations and hence of summary statistics are simulated and two kernel density estimates derived, both to be applied to the observed summary. The sequences underlying the Kiefer-Wolfowitz algorithm are taken from (the excellent optimisation book of) Spall (2003). Along with on-the-go adaptation and convergence test.

The theoretical difficulty in this extension is however that the kernel density estimator is not unbiased and thus that, rigorously speaking, the validation of the Kiefer-Wolfowitz algorithm does not apply here. On the practical side, the need for multiple starting points and multiple simulations of pseudo-samples may induce considerable time overload. Especially if  bootstrap is used to evaluate the precision of the MLE approximation. Besides normal and M/G/1 queue examples, the authors illustrate the approach on a population genetic dataset of Borneo and Sumatra orang-utans. With 5 parameters and 28 summary statistics. Which thus means using a kernel density estimator in dimension 28, a rather perilous adventure..!

future of computational statistics

Posted in Books, pictures, R, Statistics, University life with tags , , , , , , , , , , , , , , on September 29, 2014 by xi'an

I am currently preparing a survey paper on the present state of computational statistics, reflecting on the massive evolution of the field since my early Monte Carlo simulations on an Apple //e, which would take a few days to return a curve of approximate expected squared error losses… It seems to me that MCMC is attracting more attention nowadays than in the past decade, both because of methodological advances linked with better theoretical tools, as for instance in the handling of stochastic processes, and because of new forays in accelerated computing via parallel and cloud computing, The breadth and quality of talks at MCMski IV is testimony to this. A second trend that is not unrelated to the first one is the development of new and the rehabilitation of older techniques to handle complex models by approximations, witness ABC, Expectation-Propagation, variational Bayes, &tc. With a corollary being an healthy questioning of the models themselves. As illustrated for instance in Chris Holmes’ talk last week. While those simplifications are inevitable when faced with hardly imaginable levels of complexity, I still remain confident about the “inevitability” of turning statistics into an “optimize+penalize” tunnel vision…  A third characteristic is the emergence of new languages and meta-languages intended to handle complexity both of problems and of solutions towards a wider audience of users. STAN obviously comes to mind. And JAGS. But it may be that another scale of language is now required…

If you have any suggestion of novel directions in computational statistics or instead of dead ends, I would be most interested in hearing them! So please do comment or send emails to my gmail address bayesianstatistics