# Find multinominal maximum with NMaximize?

Posted 4 months ago
740 Views
|
10 Replies
|
8 Total Likes
|
 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}}
10 Replies
Sort By:
Posted 4 months 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^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 4 months 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 4 months 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*(-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 4 months 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^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 4 months 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 4 months 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[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 4 months 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.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 4 months ago
 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.