Quasi Probability method

clear all; clc;
theta = pi/4;
p = 0.3;
% initial state
psi = 1/sqrt(2)*[1;exp(1j*theta)];
rho = psi*psi';
sx = [0 1;1 0]; sy = [0 -1j;1j 0]; sz = [1 0;0 -1];
% niosy rho
nrho = (1 - 0.75*p)*rho + p/4*(sx*rho*sx' + sy*rho*sy' + sz*rho*sz');
% measurement basis
up = 1/sqrt(2)*[1;1]; down = 1/sqrt(2)*[1;-1];
n = 3000;
res1 = ones(n,1);
res2 = ones(n,1);
% i stands for times of measurement
for i = 1:n
    sum = 0;
    for j = 1:i
        num = rand;
        if num < up'*rho*up
            sum = sum + 1;
        else 
            sum = sum - 1;
        end
    end
    sum = sum / i;
    res1(i) = sum;
end
h1 = histogram(res1,'Normalization','probability'); hold on;

% Showing that the reverse only lack a constant factor, then we can have
% the same precision.
% Numerical is the same as theoretical method.
N = ((p+2)^2/(4*(1-p)^2)-cos(theta)^2)/sin(theta)^2*n;
for i = 1:N
    sum = 0;
    for j = 1:i
        flag = 0;
        measure_rand = rand;
        if measure_rand < p/(2*p + 4)
            measure_rho = sx*nrho*sx';
            flag = 1;
        elseif measure_rand < p/(p + 2)
            measure_rho = sy*nrho*sy';
            flag = 2;
        elseif measure_rand < 3*p/(2*p + 4)
            measure_rho = sz*nrho*sz';
            flag = 3;
        else
            measure_rho = nrho;
        end
        num = rand;
        if num < up'*measure_rho*up
            if flag == 0
                sum = sum + 1;
            else
                sum = sum - 1;
            end
        else 
            if flag == 0
                sum = sum - 1;
            else
                sum = sum + 1;
            end
        end
    end
    sum = (p+2)/(2-2*p)*sum / i;    
%     scatter(sum,i,25,'filled'); hold on;
    res2(i) = sum;
end
h2 = histogram(res2,'Normalization','probability');
h1.BinWidth = 0.01;
h2.BinWidth = 0.01;
legend('no error','quasi-probability');