# Why can't I get the maximum value?

   (* structural entropy for N=3 *)
S[p1_, p2_, p3_] := (
p1 + p2 + p3 == 1;
S1 = -(p1*Log[p1] + p2*Log[p2] + p3*Log[p3]);
S2 = -Log[(p1^2 + p2^2 + p3^2)];
Sstr = S1 - S2   )

NMaximize[{S[p1, p2, p3], 0 <= p1 <= 1, 0 <= p2 <= 1, 0 <= p3 <= 1, p1 + p2 + p3 == 1},
{p1, p2, p3}]


I get

{0.137934, {p1 -> 0.0320363, p2 -> 0.0223545, p3 -> 0.945609}}


while

S[0.806, 0.097, 0.097]


you can get

0.223654


So I believe 0.137934 is not the maximum value

• p1, p2, p3 are not linearly independent . you may replace p3 by 1 - p1 - p2 . Further, the value of S can be complex for 0 < {p1, p2} < 1. Therefore, there is no maximum . You must restrict p1 and p2 further to ensure a real value of S. Aug 28, 2022 at 9:34
• I can't think out other restriction for p1 and p2 Aug 28, 2022 at 9:43

Method -> "RandomSearch" seems work.

Clear[S1, S2, Sstr];
S1 = -(p1*Log[p1] + p2*Log[p2] + p3*Log[p3]);
S2 = -Log[p1^2 + p2^2 + p3^2];
Sstr = S1 - S2
NMaximize[{Sstr, 0 < {p1, p2, p3} < 1, p1 + p2 + p3 == 1}, {p1, p2,
p3}, Method -> "RandomSearch"]


The same as

Clear[S1, S2, diff, sol];
S1 = -(p1*Log[p1] + p2*Log[p2] + p3*Log[p3]);
S2 = -Log[p1^2 + p2^2 + p3^2];
diff = S1 - S2 /. p3 -> 1 - p1 - p2;
sol = NMaximize[{diff, 0 < {p1, p2} < 1, p1 + p2 < 1}, {p1, p2},
Method -> "RandomSearch"]

• Great! How did you know that? Aug 28, 2022 at 9:44
• @KarryMa Read the help document ,there are some Method in options. Aug 28, 2022 at 9:56

First have a look at what is going on and then give good starting intervals for NMaximize.

Edited There are three equal maxima.

S[p1_, p2_, p3_] = (p1 + p2 + p3 == 1;
S1 = -(p1*Log[p1] + p2*Log[p2] + p3*Log[p3]);
S2 = -Log[(p1^2 + p2^2 + p3^2)];
Sstr = S1 - S2)

Plot3D[S[p1, p2, 1 - p1 - p2], {p1, 0, 1}, {p2, 0, 1},
PlotRange -> {0, 1/2}]

Reduce[S[p1, p2, 1 - p1 - p2] \[Element] Reals, p2, Reals]

(*   0 < p1 < 1 && 0 < p2 < 1 - p1   *)

nmax = NMaximize[{S[p1, p2, 1 - p1 - p2],
0 < p2 < 1 - p1}, {{p1, 0, .1}, {p2, 0, .1}}]

(*   {0.223655, {p1 -> 0.096693, p2 -> 0.096693}}   *)

nmax = NMaximize[{S[p1, p2, 1 - p1 - p2],
0 < p2 < 1 - p1}, {{p1, .8, 1}, {p2, 0, 1}}]

(*   {0.223655, {p1 -> 0.806614, p2 -> 0.0966931}}   *)

nmax = NMaximize[{S[p1, p2, 1 - p1 - p2],
0 < p2 < 1 - p1}, {{p1, 0, 1}, {p2, .8, 1}}]

(*   {0.223655, {p1 -> 0.096693, p2 -> 0.806614}}   *)

• It seems that it is not the answer Aug 28, 2022 at 9:49
• Corrected my answer. Aug 28, 2022 at 18:49