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?

Improve this question
$\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
    Commented Sep 3, 2020 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
    Commented Sep 4, 2020 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
    Commented Sep 4, 2020 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
    Commented Sep 4, 2020 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
    Commented Sep 4, 2020 at 4:31

1 Answer 1

Reset to default
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 *)
Improve this answer
$\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
    Commented Sep 4, 2020 at 4:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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