Group Abstract

Message Boards

WOLFRAM COMMUNITY

8.2K Views

10 Replies

8 Total Likes

View groups...

Follow this post

Share this post:

GROUPS:

Mathematics Wolfram Language Optimization Numerical Computation

Find multinominal maximum with NMaximize?

Alex Graham

Posted 7 years ago

POSTED BY: Alex Graham

10 Replies

Sort By:

Frank Kampas

Frank Kampas, Physicist at Large Consulting

Posted 7 years ago

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

Alex Graham

Posted 7 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

Frank Kampas

Frank Kampas, Physicist at Large Consulting

Posted 7 years ago

Turns out only one of the 50 roots is nonzero. In[35]:= res2a = Reduce[{obj == -n!/((94)!4!)p0^(94)ap10^4p1*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.

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

Frank Kampas

Frank Kampas, Physicist at Large Consulting

Posted 7 years ago

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.145110^-6, \[Mu][2.] -> -2.2511910^-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}

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

Frank Kampas

Frank Kampas, Physicist at Large Consulting

Posted 7 years ago

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[musineqcons == 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^4p1 p2, {p0 + p1 + p100 + p2 + ap10 == 1, (1/100)Y(-10ap10 + 30p1 + 40p2 + 100p100) == 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^4p1*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

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

Alex Graham

Posted 7 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^4p1p2, p0 + p1 + p100 + p2 + ap10 == 1 && p0 >= 0 && p1 >= 0 && p2 >= 0 && ap10 >= 0 && p100 >= 0 && (1/100)Y(-10ap10 + 30p1 + 40p2 + 100p100) == sto}, {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.3947710^-9, ap10 -> 0.0437649}}} If you dont mind please take a look at a bit more complicated example. Attachments:* Maximums.nb

POSTED BY: Alex Graham

Frank Kampas

Frank Kampas, Physicist at Large Consulting

Posted 7 years ago

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^4p1* p2, {p0 + p1 + p100 + p2 + ap10 == 1, (1/100)Y(-10ap10 + 30p1 + 40p2 + 100p100) == 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.}}

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

Frank Kampas

Frank Kampas, Physicist at Large Consulting

Posted 7 years ago

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(-10ap10 + 30p1 + 40p2 + 100p100) == 1000}, {p0, p1, p2, p100, ap10}]; In[16]:= Quiet[ NMaximize[{n!/((94)!4!)p0^(94)ap10^4p1p2 , #}, {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.8752410^-10, {p0 -> 0.997462, p1 -> 0.0011485, p2 -> 0.000577472, p100 -> 0.000459685, ap10 -> 0.000352236}} \|\| {3.7625610^-8, {p0 -> 0.996667, p1 -> 0.000555751, p2 -> 0.000648745, p100 -> 0.000715146, ap10 -> 0.00141369}} \|\| {2.8726910^-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}}

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

Alex Graham

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^4p1p2, p0 + p1 + p100 + p2 + ap10 == 1 && p0 >= 0 && p1 >= 0 && p2 >= 0 && ap10 >= 0 && p100 >= 0 && (1/100)Y(-10ap10 + 30p1 + 40p2 + 100p100) == sto}, {p0, p1, p2, p100, ap10}, Method -> "DifferentialEvolution"], {o, o1, o2}] Transpose[ Join[{Table[sto, {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`"}

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

Frank Kampas

Frank Kampas, Physicist at Large Consulting

Posted 7 years ago

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^4p1p2, p0 + p1 + p100 + p2 + ap10 == 1 && p0 >= 0 && p1 >= 0 && p2 >= 0 && ap10 >= 0 && p100 >= 0 && (1/100)Y(-10ap10 + 30p1 + 40p2 + 100p100) == 1000}, {p0, p1, p2, p100, ap10}, Method -> "DifferentialEvolution"] Out[33]= {0.0251849, {p0 -> 0.940801, p1 -> 0.00797022, p2 -> 0.00746355, p100 -> 2.6215710^-9, ap10 -> 0.0437649}}

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

Reply Preview

Attachments

Remove Add a file to this post

Follow this discussion

or Discard

Feedback