Message Boards Message Boards

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

Find multinominal maximum with NMaximize?

Posted 6 years ago

Hello guys. This is the table with the result that should be got https://www.dropbox.com/s/rohvusatedkfh3r/Screenshot_20180725-184452.png?dl=0 The mysticism thing is that it found a part of correct values only with NMaximize. With 3000-26000 everithing is OK but with others a strange thing is happening. For example, if I take 3000 and more I get the right answer

Clear["Global`*"]
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) == 
     3000}, {p0, p1, p2, p100, ap10}]
{0.0272692, {p0 -> 0.94, p1 -> 0.00999995, p2 -> 0.00999991, 
  p100 -> 5.71366*10^-8, ap10 -> 0.0400001}}

So I took all values from -6000 to 29000 and get this. Most of the values are correct, but some of the definatele wrong.

{
RowBox[{"-", "6000"}], "0.00273003019908044`"},
{
RowBox[{"-", "5000"}], "4.002433792920723`*^-12"},
{
RowBox[{"-", "4000"}], "9.420778019991495`*^-80"},
{
RowBox[{"-", "3000"}], "0.010805254786454756`"},
{
RowBox[{"-", "2000"}], "0.014835126377371699`"},
{
RowBox[{"-", "1000"}], "1.1255118610736836`"},
{"0", "1.0083957979097389`"},
{"1000", "3.3664295470533197`*^-22"},
{"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`"},
{"11000", "0.01467181808621124`"},
{"12000", "0.013260174294049413`"},
{"13000", "0.011983123654999555`"},
{"14000", "0.01082795065726872`"},
{"15000", "0.00978312998048758`"},
{"16000", "6.311707230722985`*^-18"},
{"17000", "0.007983743159662302`"},
{"18000", "0.007211133780524158`"},
{"19000", "0.00651261660134752`"},
{"20000", "0.005881151176722276`"},
{"21000", "0.005310359883946591`"},
{"22000", "0.004794466105647533`"},
{"23000", "0.004328238122009399`"},
{"24000", "0.003906938191582962`"},
{"25000", "0.0035262763453895838`"},
{"26000", "0.0031823684622467223`"},
{"27000", "4.580651706257595`*^-15"},
{"28000", "2.32433808006728`*^-15"},
{"29000", "1.5695861700229384`*^-14"}

That is more interesting with reduce function gives different results.

Clear["Global`*"]
n = 100; 
Y = 1000000; 
r = Reduce[
   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) == 
     3000, {ap10, p0, p1, p2, p100}, Reals, Backsubstitution -> True];
q = NMaximize[{n!/((94)!*4!)*p0^(94)*ap10^4*p1*p2, r}, {p0, p1, p2, 
   p100, ap10}]

{0.0230288, {p0 -> 0.940867, p1 -> 0.00763999, p2 -> 0.00891332, p100 -> 0.00127333, ap10 -> 0.0413065}} FInd Maximum also dosent work

Clear["Global`*"]
Y = 1000000;
n = 100; 
q = FindMaximum[{n!/((94)!*4!)*p0^(94)*p1*p2*p10^4, 
   p0 + p1 + p100 + p2 + p10 == 1 && p0 >= 0 && p1 >= 0 && p2 >= 0 && 
    p100 >= 0 && 
    p10 >= 0 && (1/100)*Y*(-10*p10 + 30*p1 + 40*p2  + 100*p100) == 
     3000}, {{p0, 0.91}, {p1, 0.01387}, {p2, 0.016}, {p100, 
    0.021}, {p10, 0.35}}]
{4.61652*10^-31, {p0 -> 0.412788, p1 -> 0.0133285, p2 -> 0.0185802, 
  p100 -> 0.0428179, p10 -> 0.512485}}

So what is that, wrong table at dropbox or I need to use a different code? Multiple searches also dont do a good thing,http://community.wolfram.com/groups/-/m/t/1164680

n = 100; 
Y = 1000000; 
iMin[-n!/((94)!*4!)*p0^(94)*p1*p2*p10^4, 
 List @@ (p0 + p1 + p2 + p10 + p100 == 1 &&
    p0 >= 0 && p1 >= 0 && p2 >= 0 && p100 >= 0 && 
    p10 >= 0 && (1/100)*Y*(-10*p10 + 30*p1 + 40*p2  + 100*p100) == 
     3000),
 Thread[{{p0, p1, p2, p100, p10}, 0, 1}], 10, 0]
{-1.57794*10^-44, {p0 -> 0.289568, p1 -> 0.0329513, p2 -> 0.0407228, 
  p100 -> 0.0368193, p10 -> 0.599939}}
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 6 years ago

Thank you for your reply. This solution dosent fits me as I need to find around 100 maximums. And in a more complicated example it has much more solutions and it takes too much time. What really helps is to add extra conditions. We know that the probability cannot be as electron diameter, and I add extra condition 1>=n!/((94)!4!)p0^(94)ap10^4p1*p2>=1/10^6 and use reduce function not for all condition. What is more interesting, multiple random seraches gives the same bad result even with this extra condition. I have some question for you. I'm going to use multiple random serches and kurh tucker equations in my research. You developed these methods, and I will make a reference link of it. But maybe you publish these methodes someware else? In this case, it's better to give a link at some of your article directly.

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 6 years ago

Thank you for your answer. I've made some experiments and found that it gives good results with adding scalingfactor->1. But if equations becomes more complicated it gives strange results again. So what does scalefactor do and which paramenter should I use or maybe I should add some more methodes? According to iMin with global it gives better results giving the zero but the right answer is 0.025. It would be amazing if NMaximize function will be the single method of finding this naximum. But even NMaximize cant cope with more complicated equations((

n = 100; 
Y = 1000000; 
st = 1000;
o1 = 1;
o2 = 1;

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", "ScalingFactor" -> 1}], {o, 
    o1, o2}];
{{0.0251849, {p0 -> 0.940801, p1 -> 0.00797022, p2 -> 0.00746355, 
   p100 -> 4.39477*10^-9, ap10 -> 0.0437649}}}

If you dont mind please take a look at a bit more complicated example.

Attachments:
POSTED BY: Alex Graham

using some experimental code I've written which uses ParametricIPOPTMinimize to do minimization (and changing the sign of the objective function). I got an objective function of 0 for the case where the second constraint has a right-hand side of 1000.

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

In[24]:= iMin[global, -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}, 
 Thread[{{p0, p1, p2, p100, ap10}, 0, 1}], 1, 0]

During evaluation of In[24]:= solving time = 0.00369322

During evaluation of In[24]:= {{Solve_Succeeded,1}}

Out[24]= {{0., {0.179392, 0.060826, 0.018392, 0.045031, 0.696359}, 
  0.}}
POSTED BY: Frank Kampas

You might get better results by using LogicalExpand on the result of Reduce on the constraints and taking the best result.

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

In[15]:= r1 = 
  Reduce[{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}];

In[16]:= Quiet[
 NMaximize[{n!/((94)!*4!)*p0^(94)*ap10^4*p1*p2 , #}, {p0, p1, p2, 
     p100, ap10}, Method -> "DifferentialEvolution"] & /@ 
  LogicalExpand[r1]]

Out[16]= {0., {p0 -> 0.996973, p1 -> 0.00289558, p2 -> 0, 
   p100 -> 0.000131325, ap10 -> 0}} || {0., {p0 -> 0.998052, 
   p1 -> 0.00135491, p2 -> 0, p100 -> 0.000593527, 
   ap10 -> 0}} || {0., {p0 -> 0.630183, p1 -> 0.0949542, p2 -> 0, 
   p100 -> 0, ap10 -> 0.274863}} || {2.87524*10^-10, {p0 -> 0.997462, 
   p1 -> 0.0011485, p2 -> 0.000577472, p100 -> 0.000459685, 
   ap10 -> 0.000352236}} || {3.76256*10^-8, {p0 -> 0.996667, 
   p1 -> 0.000555751, p2 -> 0.000648745, p100 -> 0.000715146, 
   ap10 -> 0.00141369}} || {2.87269*10^-9, {p0 -> 0.9975, 
   p1 -> 0.000357294, p2 -> 0.000416864, p100 -> 0.000816955, 
   ap10 -> 0.000908887}} || {0.0251849, {p0 -> 0.940801, 
   p1 -> 0.00797022, p2 -> 0.00746356, p100 -> 1.02811*10^-12, 
   ap10 -> 0.0437649}}
POSTED BY: Frank Kampas
Posted 6 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

Group Abstract Group Abstract