After: p = p ./ numel(I), sum(p) is not 1.0, but 0.333. You do not want to divide by numel(I), but by numel(Red), or equivalently: size(I,1)*size(I,2).
I = rand(640, 480, 3);
Red = I(:,:,1);
Green = I(:,:,2);
Blue = I(:,:,3);
p = imhist(Red); % Histogram
p(p == 0) = []; % remove zero entries in p
p = p ./ numel(Red); % normalize p so that sum(p) is 1,0
Er = round(-sum(p .* log2(p)), 3)
p = imhist(Blue); % Histogram
p(p == 0) = []; % remove zero entries in p
p = p ./ numel(Blue); % normalize p so that sum(p) is one.
Eb = round(-sum(p .* log2(p)), 3)
You see, that the results are not 0 for random inputs.
If you want to collect the data, follow Scott's advice:
nFile = numel(myFiles);
Er = zeros(1, nFile);
Eb = zeros(1, nFile);
for k = 1 : nFile
...
Er(k) = round(-sum(p .* log2(p)), 3);
...
Eb(k) = round(-sum(p .* log2(p)), 3);
end
percentage = sum(Er > Eb) / nFile * 100;
fprintf('Percentage of images with red entropy > blue entropy: %f\n', percentage);