0
$\begingroup$

I'm trying to use NMaximize to show that entropy is maximized when the distribution is uniform. It seems to work when the random variable is valued, but breaks at 4+ variables:

p = {p1, p2, p3, p4};
c = Map[0 <= # <= 1 &, p];
NMaximize[{Sum[-Log[p[[i]]]*p[[i]], {i, 1, Length[p]}], Total[p] == 1,
   Splice[c]}, p]

If p is changed to {p1, p2, p3}, then the result is p1 = p2 = p3 = 1/3 as expected. But for 4 probabilities, I get the following:

{1.30778, {p1 -> 0.3323, p2 -> 0.269063, p3 -> 0.0975481, 
  p4 -> 0.301089}}

Which seems clearly wrong. Am I using NMaximize incorrectly?

$\endgroup$
8
  • 2
    $\begingroup$ Weirdly this works using the $L_1$ norm i.e Abs[p1] + Abs[p2] + Abs[p3] + Abs[p4] by doing - NMaximize[{Total[-# Log[#] & /@ p], Splice[c], Norm[p, 1] == 1}, p] and gives the result: {1.38629, {p1 -> 0.250001, p2 -> 0.25, p3 -> 0.25, p4 -> 0.25}} but note the slight error of .000001 which violates the constraint. $\endgroup$
    – flinty
    Sep 3 '20 at 23:56
  • 2
    $\begingroup$ If you use the "RandomSearch" method then it gives the right result straight away: NMaximize[{Total[-# Log[#] & /@ p], Splice[c], Total[p] == 1}, p, Method -> "RandomSearch"] $\endgroup$
    – flinty
    Sep 4 '20 at 0:02
  • $\begingroup$ Just tell NMaximize about the symmetry NMaximize[{Sum[-Log[p[[i]]]*p[[i]], {i, 1, Length[p]}], Total[p] == 1, p1 == p2 == p3 == p4, c} // Flatten, p] . This works also with Maximize Maximize[{Sum[-Log[p[[i]]]*p[[i]], {i, 1, Length[p]}], Total[p] == 1, p1 == p2 == p3 == p4, c} // Flatten, yields {Log[4], {p1 -> 1/4, p2 -> 1/4, p3 -> 1/4, p4 -> 1/4}} $\endgroup$
    – Akku14
    Sep 4 '20 at 4:25
  • $\begingroup$ @flinty The L1 norm trick doesn't seem to generalize if I add more parameters, but RandomSearch seems to work well consistently. $\endgroup$
    – rgrinberg
    Sep 4 '20 at 4:31
  • $\begingroup$ @Akku14 the symmetry is something I knew about the nature of the solution in this case. It's not going to help me in other situations where the probabilities are going to be different $\endgroup$
    – rgrinberg
    Sep 4 '20 at 4:31
1
$\begingroup$
Clear["Global`*"]

p = {p1, p2, p3, 1 - p1 - p2 - p3};

sol = Solve[Thread[D[Total[-p*Log[p]], {Most@p}] == 0], Most@p]

(* {{p1 -> 1/4, p2 -> 1/4, p3 -> 1/4}} *)

p /. sol[[1]]

(* {1/4, 1/4, 1/4, 1/4} *)

Total[-p*Log[p]] /. sol[[1]]

(* Log[4] *)

EDIT: For larger number of probabilities

n = 50;

p = pr /@ Range[n - 1];
p = Append[p, 1 - Total[p]];

sol = Solve[Thread[D[Total[-p*Log[p]], {Most@p}] == 0], Most@p];

(p /. sol[[1]] // Union) == {1/n}

(* True *)

(Total[-p*Log[p]] /. sol[[1]]) == Log[n]

(* True *)
$\endgroup$
1
  • 1
    $\begingroup$ That works alright, but it's a little manual. I would prefer to use mathematica's native optimization functions. $\endgroup$
    – rgrinberg
    Sep 4 '20 at 4:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.