Q:

NMaximize seems to be returning the wrong result

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?

A:

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 *)

