Why does my simulation of the binary symmetric channel disagree with theory?

Question

Could someone please answer me why the simulated performance (using matlab) of a Hamming code over a binary symmetric channel does not fit exacly with the theoretical? Here is my code:

clear all

n=15;
k=11;
I=eye(k);

%define the parity matrix for the hamming code 
P=[1 1 0 0;1 0 1 0;0 1 1 0;1 0 0 1;0 1 0 1;0 0 1 1;1 1 1 0;1 1 0 1;1 0 1 1;0 1 1 1;1 1 1 1];

%The generator matrix combination of parity and identity matrix
G=[I,P];

%Define parity check matrix
HI=eye(n-k);
H=[P',HI];

%verify that rows of H represent the dual space for the rows of G
%(verification = 0)
verification=rem(G*H',2);

UoB=1043162;

%convert UoB number to binary
ci=dec2bin(UoB);
ci=ci-('0');

%take the 11 rightmost digits
cw=ci(:,10:20);
cw=rem(cw*G,2);

%monte carlo simulation index
mc=10000;
index=1;
Pe_matrix=zeros(21,mc);
Pe_vector=zeros(1,21);

In the procedure demonstrated bellow random error patterns error=rand(1,15)<p; are generated for 21 different crossover probabilities.Consequently these patterns are subtracted from the original codeword c=xor(cw,error); in order to find the errored codeword (in any place where the error pattern has an 1 the codeword will flip its bit. After that, the probability of error for the current error patern is calculated Pe=(sum(xor(cw,c)))./k; by summing up the ones returned (the places where the errored codeword and the original differ) from the subtraction of the errored codeword and the original one, and by dividing them by the total length of the codeword k=numel(c);.The Pe_matrix stores all the Pe values for each of the cross over probabilities and each of the mc iterations.

for i=1:mc

    for p=0:0.01:.2   
        error=rand(1,15)<p;
        c=xor(cw,error);
        k=numel(c);
        Pe=(sum(xor(cw,c)))./k;
        Pe_matrix(index,i)=Pe;
        index=index+1;
    end
    index=1;

end

In order to obtain smoother results Monte Carlo simulations were performed. In each of the Monte Carlo iterations,each row of the Pe_matrix has the Pe values for the different crossover probabilities.Now by summing up the elements of each row and by dividing them with the total number of the Monte Carlo simulations we acquire a strong estimate of the BER for the Binary symmetric channel.

for j=1:21
    Pe_vector(1,j)=sum(Pe_matrix(j,:))/mc;
end

%In the code demonstrated bellow the theoretical performance of the channel
%is calculated in order to compare it with the simulated.

for pe=0:0.01:.2
    index=index+1;
    Pblkerror=1-(1-pe)^15-15*pe*(1-pe)^14;
    Pd(1,index)=(3/15)*Pblkerror;
end

%finally by creating a legend one can compare the theoretical with the
%simulated performance by the depicted BER graphs.
semilogy(Pe_vector,'r-o')
hold on;
semilogy(Pd,'k-*')
hold on;
legend('simulation','theoretical performance');
ylabel('Probability of error,Pe');
xlabel('cross over probability');
title('BsC simulation');
grid on;