share/cobyla/ex/die.mac

   1 /*
   2  * Maximum entropy estimate.  Based on an example from Kostas
   3  * Oikonomou.
   4  *
   5  * Given a possibly unfair die where each side n has a probability
   6  * of occurrence of xn, find the most likely value if it is known
   7  * that the average value is av.
   8  *
   9  * Let the probabilities be x1, x2, x3, x4, x5, x6.  The criterion for
  10  * most likely is maximizing the entropy H(X) = sum(-x*log(x), all X).
  11  * We are subject to the constraints x1+x2+x3+x4+x5+x6 = 1, of course,
  12  * and to x1+2*x2+3*x3+4*x4+5*x5+6*x6 = av, since we know the average
  13  * value.  And, of course, since these are probabilities, 0 <= xn <= 1.
  14  *
  15  * This example shows how to get equality constraints when fmin_coblya
  16  * normally wants inequality constraints.  If we want g(x) = 0, we can use
  17  * the inequality constraints g(x) >= 0 and -g(x) >= 0.
  18  *
  19  */
  20
  21 HH(X) := -xreduce("+",map(lambda([x],x*log(x)),X));
  22
  23 /*
  24  * Although we want xn >= 0, the algorithm sometimes sets xn = 0, which
  25  * causes problems when we evaluate xn*log(xn).  So, instead of xn >= 0,
  26  * we use xn >= 1d-9.  (That is, xn - 1d-9 >= 0.)
  27  */
  28 die(av) :=
  29   fmin_cobyla(-HH([x1,x2,x3,x4,x5,x6]),
  30               [x1,x2,x3,x4,x5,x6],
  31               makelist(1/6,k,1,6),
  32               iprint=1,
  33               constraints= [
  34                 x1 >= 1d-9, x2 >= 1d-9, x3 >= 1d-9, x4 >= 1d-9, x5 >= 1d-9, x6 >= 1d-9,
  35                 x1<=1, x2<=1, x3<=1, x4<=1, x5<=1, x6<=1,
  36                 x1+x2+x3+x4+x5+x6-1 = 0,
  37                 1*x1+2*x2+3*x3+4*x4+5*x5+6*x6=av]);
  38
  39
  40 /*
  41  * If the faces are all equally likely, the average would be 7/2.
  42  * Here's what the maximum-entropy result would be:
  43  */
  44
  45 die(7/2);
  46 /*
  47    Normal return from subroutine COBYLA
  48
  49    NFVALS =  146   F =-1.791759E+00    MAXCV = 1.110223E-15
  50    X = 1.666669E-01   1.666665E-01   1.666660E-01   1.666672E-01   1.666669E-01
  51        1.666665E-01
  52
  53 [[x1 = .1666669204062052,x2 = .1666665386187796,
  54           x3 = .1666659607581757,x4 = 0.166667233111882,
  55           x5 = 0.166666894995565,x6 = .1666664521093927],-1.791759469225061,
  56          146, 0]
  57
  58 We see that the probabilities are close to 1/6, as expected.
  59 */
  60
  61 /*
  62  * A more interesting test is if the average is 4.5
  63  */
  64 die(4.5);
  65 /*
  66    Normal return from subroutine COBYLA
  67
  68    NFVALS =  160   F =-1.613581E+00    MAXCV = 1.154632E-14
  69    X = 5.435361E-02   7.877099E-02   1.141600E-01   1.654467E-01   2.397745E-01
  70        3.474942E-01Evaluation took 0.0500 seconds (0.0500 elapsed) using 1.443 MB.
  71 [[x1 = .05435294488099255,x2 = .07877230426143854,
  72           x3 = 0.1141598558557,x4 = .1654458115488451,x5 = .2397748678844901,
  73           x6 = .3474942155685333],-1.613581098146268,
  74          142, 0]
  75 */
  76
  77 /*
  78  * The obvious values for die(1) is x1=1 all others are 0.
  79  * For die(6), we must have x6 = 1, and all other are 0.
  80  */
  81
  82 die(1);
  83 /*
  84
  85    Normal return from subroutine COBYLA
  86
  87    NFVALS =   58   F =-1.491961E-08    MAXCV = 8.823530E-10
  88    X = 1.000000E+00   1.176470E-10   1.176471E-10   1.176471E-10   1.176471E-10
  89        1.176471E-10Evaluation took 0.1100 seconds (0.1400 elapsed) using 4.477 MB.
  90 [[x1 = .9999999985294118,x2 = 1.176470448838174e-10,
  91           x3 = 1.176470617973713e-10,x4 = 1.176470548584774e-10,
  92           x5 = 1.176470500879878e-10,x6 = 1.176470604963287e-10],
  93          -1.4919606548618195e-8,
  94          58, 0]
  95 */
  96
  97 die(6);
  98 /*
  99
 100    Normal return from subroutine COBYLA
 101
 102    NFVALS =   56   F =-3.569971E-08    MAXCV = 6.818184E-10
 103    X = 3.181816E-10   3.181817E-10   3.181817E-10   3.181817E-10   3.181816E-10
 104        1.000000E+00
 105 [[x1 = 3.18181639022419e-10,x2 = 3.181816702474416e-10,
 106           x3 = 3.181817014724642e-10,x4 = 3.181816546349303e-10,
 107           x5 = 3.181815921848852e-10,x6 = .9999999990909101],
 108          -3.5699705929961904e-8,
 109          56, 0]
 110 */