Group Abstract Group Abstract

Message Boards Message Boards

0
|
8.4K Views
|
10 Replies
|
8 Total Likes
View groups...
Share
Share this post:

Find multinominal maximum with NMaximize?

Posted 7 years ago
POSTED BY: Alex Graham
10 Replies

I wrote up solving the Kuhn-Tucker equations with Reduce some years ago in the Mathematica Journal.

http://www.mathematica-journal.com/issue/v9i4/contents/Tricks9-4/Tricks9-4_2.html

There are some typos in the article but you can use it for a reference.

POSTED BY: Frank Kampas
Posted 7 years ago
POSTED BY: Alex Graham

Turns out only one of the 50 roots is nonzero.

In[35]:= res2a = 
 Reduce[{obj == -n!/((94)!*4!)*p0^(94)*ap10^4*p1*p2, res1}, {obj, p0, 
     p1, p2, p100, ap10}, Backsubstitution -> True][[All, 8]] // N

Out[35]= obj == 0. || obj == 0. || obj == 0. || obj == 0. || 
 obj == 0. || obj == 0. || obj == 0. || obj == 0. || obj == 0. || 
 obj == 0. || obj == 0. || obj == 0. || obj == 0. || obj == 0. || 
 obj == 0. || obj == 0. || obj == 0. || obj == 0. || obj == 0. || 
 obj == 0. || obj == 0. || obj == -0.0251849 || obj == 0. || 
 obj == 0. || obj == 0. || obj == 0. || obj == 0. || obj == 0. || 
 obj == 0. || obj == 0. || obj == 0. || obj == 0. || obj == 0. || 
 obj == 0. || obj == 0. || obj == 0. || obj == 0. || obj == 0. || 
 obj == 0. || obj == 0. || obj == 0. || obj == 0. || obj == 0. || 
 obj == 0. || obj == 0. || obj == 0. || obj == 0. || obj == 0. || 
 obj == 0. || obj == 0. || obj == 0.
POSTED BY: Frank Kampas

Looking at the other solutions, the 22nd solution has an objective value(as a maximum problem) of 0.0251849

In[28]:= res2rl = List @@ (ToRules /@ res2);

In[29]:= Head /@ res2rl

Out[29]= {List, List, List, List, List, List, List, List, List, List, \
List, List, List, List, List, List, List, List, List, List, List, \
List, ToRules, ToRules, ToRules, ToRules, ToRules, ToRules, ToRules, \
ToRules, ToRules, ToRules, ToRules, ToRules, ToRules, ToRules, \
ToRules, ToRules, ToRules, ToRules, ToRules, ToRules, ToRules, \
ToRules, ToRules, ToRules, ToRules, ToRules, ToRules, ToRules, ToRules}

In[30]:= obj /. res2rl[[Range[22]]] // N

Out[30]= {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., \
0., 0., 0., 0., 0., 0., -0.0251849}

In[31]:= res2rl[[22]] // N

Out[31]= {\[Mu][1.] -> 0., \[Mu][3.] -> 0., \[Mu][5.] -> 
  0., \[Lambda][1.] -> 2.51635, \[Lambda][2.] -> 
  2.1451*10^-6, \[Mu][2.] -> -2.25119*10^-10, \[Mu][4.] -> 4.66145, 
 obj -> -0.0251849, p0 -> 0.940801, p1 -> 0.00797022, 
 p2 -> 0.00746356, p100 -> 3.57333*10^-17, ap10 -> 0.0437649}
POSTED BY: Frank Kampas

If I symbolically solve the Kuhn-Tucker equations for the problem with the second constraint = 1000, I get 51 solutions, the first of which has an objective function value of 0

In[1]:= KTEqs[obj_, cons_List, vars_] :=
 Module[{consconvrule = {GreaterEqual[x_, y_] -> LessEqual[y - x, 0], 
     Equal[x_, y_] -> Equal[x - y, 0], 
     LessEqual[x_, y_] -> LessEqual[x - y, 0],
     LessEqual[lb_, x_, ub_] -> LessEqual[(x - lb) (x - ub), 0], 
     GreaterEqual[ub_, x_, lb_] -> LessEqual[(x - lb) (x - ub), 0]} 
    , stdcons, eqcons, ineqcons, lambdas, mus, lagrangian, eqs1, eqs2,
    eqs3, alleqns, allvars },
  stdcons = cons /. consconvrule;
  eqcons = Cases[stdcons, Equal[_, 0]][[All, 1]];
  ineqcons = Cases[stdcons, LessEqual[_, 0]][[All, 1]];
  lambdas = Array[\[Lambda], Length[eqcons]];
  mus = Array[\[Mu], Length[ineqcons]];
  lagrangian = obj + lambdas.eqcons + mus.ineqcons;
  eqs1 = Thread[ D[lagrangian, {vars}] == 0];
  eqs2 = Thread[mus >=   0];
  eqs3 = Thread[mus*ineqcons == 0];
  alleqns = Join[eqs1, eqs2, eqs3, stdcons];
  allvars = Join[vars, lambdas, mus];
  {alleqns, allvars}
  ]

In[2]:= n = 100;
Y = 1000000;

In[4]:= kts1 = 
  KTEqs[-n!/((94)!*4!)*p0^(94)*ap10^4*p1*
    p2, {p0 + p1 + p100 + p2 + ap10 == 
     1, (1/100)*Y*(-10*ap10 + 30*p1 + 40*p2 + 100*p100) == 1000, 
    p0 >= 0, p1 >= 0, p2 >= 0, p100 >= 0, ap10 >= 0}, {p0, p1, p2, 
    p100, ap10}];

In[5]:= res1 = 
  Reduce[Sequence @@ kts1, Reals, Backsubstitution -> True];

In[6]:= Length[res1]

Out[6]= 51

In[7]:= res2 = 
  Reduce[{obj == -n!/((94)!*4!)*p0^(94)*ap10^4*p1*p2, res1}, {obj, p0,
     p1, p2, p100, ap10}];

In[8]:= Length[res2]

Out[8]= 51

In[9]:= res2[[1]]

Out[9]= \[Lambda][1] == 0 && \[Lambda][2] == 0 && \[Mu][1] == 
  0 && \[Mu][2] == 0 && \[Mu][3] == 0 && \[Mu][4] == 0 && \[Mu][5] == 
  0 && obj == 0 && p0 == 0 && p1 == 101/800 && p2 == 101/1000 && 
 p100 == 0 && ap10 == 3091/4000

In[10]:= N[res2[[1]]]

Out[10]= \[Lambda][1.] == 0. && \[Lambda][2.] == 0. && \[Mu][1.] == 
  0. && \[Mu][2.] == 0. && \[Mu][3.] == 0. && \[Mu][4.] == 
  0. && \[Mu][5.] == 0. && obj == 0. && p0 == 0. && p1 == 0.12625 && 
 p2 == 0.101 && p100 == 0. && ap10 == 0.77275
POSTED BY: Frank Kampas
Posted 7 years ago
Attachments:
POSTED BY: Alex Graham
POSTED BY: Frank Kampas
POSTED BY: Frank Kampas
Posted 7 years ago

Thank you for your reply. It improves the calculations but not for all values.

n = 100; 
Y = 1000000; 
st = 1000;
o1 = -10;
o2 = 10;
q = Table[
  NMaximize[{n!/((94)!*4!)*p0^(94)*ap10^4*p1*p2, 
    p0 + p1 + p100 + p2 + ap10 == 1 && p0 >= 0 && p1 >= 0 && p2 >= 0 && 
     ap10 >= 0 && 
     p100 >=  0 && (1/100)*Y*(-10*ap10 + 30*p1 + 40*p2  + 100*p100) ==
       st*o}, {p0, p1, p2, p100, ap10}, 
   Method -> "DifferentialEvolution"], {o, o1, o2}]
Transpose[
  Join[{Table[st*o, {o, o1, o2}], q[[All, 1]]}, 1]] // MatrixForm
RowBox[{"-", "10000"}], "0.00020081432276805312`"},
{
RowBox[{"-", "9000"}], "0.0004119921224652788`"},
{
RowBox[{"-", "8000"}], "1.8936699121839722`*^-18"},
{
RowBox[{"-", "7000"}], "0.0008268121548938624`"},
{
RowBox[{"-", "6000"}], "0.00273003019908044`"},
{
RowBox[{"-", "5000"}], "4.002433792920723`*^-12"},
{
RowBox[{"-", "4000"}], "0.007309802589918361`"},
{
RowBox[{"-", "3000"}], "0.010805254786454756`"},
{
RowBox[{"-", "2000"}], "0.014835125565673924`"},
{
RowBox[{"-", "1000"}], "1.1255118610736836`"},
{"0", "1.0083957979097389`"},
{"1000", "0.025184932372196375`"},
{"2000", "1.0105517857633675`"},
{"3000", "0.027269195370073815`"},
{"4000", "0.026831313283137287`"},
{"5000", "0.025659720271071166`"},
{"6000", "0.02396928739747252`"},
{"7000", "0.021955363082878255`"},
{"8000", "0.019861954802262567`"},
{"9000", "0.017956439320638243`"},
{"10000", "0.01623208221174882`"}
POSTED BY: Alex Graham

NMaximize doesn't always give the global maximum. I often get better results with DifferentialEvolution

In[31]:= n = 100;
Y = 1000000;
q = NMaximize[{n!/((94)!*4!)*p0^(94)*ap10^4*p1*p2, 
   p0 + p1 + p100 + p2 + ap10 == 1 && p0 >= 0 && p1 >= 0 && p2 >= 0 &&
     ap10 >= 0 && 
    p100 >= 0 && (1/100)*Y*(-10*ap10 + 30*p1 + 40*p2 + 100*p100) == 
     1000}, {p0, p1, p2, p100, ap10}, 
  Method -> "DifferentialEvolution"]

Out[33]= {0.0251849, {p0 -> 0.940801, p1 -> 0.00797022, 
  p2 -> 0.00746355, p100 -> 2.62157*10^-9, ap10 -> 0.0437649}}
POSTED BY: Frank Kampas
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard