Construct i.i.d. Bernoulli variables with MATLAB

Hello,

As known from the probability laws, the following formula is true for i.i.d. variables :

On the following code, I am trying to construct 5 i.i.d. Bernoulli variables with parameter p, and verify the above formula from a large number of experiments. However, as can be seen from the printed matrices on the end of the script, the above formula fails.

I know that the reason is the lack of independence between the random variables, however I have searched a lot and I have not found any solutions. Instead, I found many answers suggesting exactly the method I use here, i.e. "diag( rand(n,1) < p )". Where am I wrong?

rng(1);
n = 4;
p = 0.1;
experiments = 100000;

for e = 1 : experiments %e different draws of every random variable, in order to obtain enough samples
   for ell=1:5  %5 different Psi variables
      Psi_matrix(:,ell,e) = ( rand(n,1) < p );
   end
end

mean1 = mean(Psi_matrix, 3);
prod1 = mean1(:,1).*mean1(:,2).*mean1(:,3).*mean1(:,4).*mean1(:,5); %product of expectations

prod0 = zeros(n,experiments);
for e=1:experiments
    prod0(:,e) = Psi_matrix(:,1,e).*Psi_matrix(:,2,e).*Psi_matrix(:,3,e).*Psi_matrix(:,4,e).*Psi_matrix(:,5,e);
end
mean0 = mean(prod0, 2);  %expectation of product

prod1  %Check the printed matrices - prod1 is unequal to mean0
mean0

[Matt J]

Why diag?

[Georgios Apostolakis]

On my implementation $\Psi^{(i)}$ were diagonal matrices with i.i.d. Bernoulli variables on their main diagon. However, that's not of significant importance and doesn't change the problem, so I edited the question in order that are random vectors. Thanks for the comment.

Answer 1

[Paul]

Following code illustrates the expected result with p = 0.4. I wasn't sure why the Question was using 3-dimensional matrices, so they are not used here.

rng(1);
p = 0.4;
experiments = 100000;
psi = rand(experiments,5) < p;
prod(mean(psi,1))

ans = 0.0102

mean(prod(psi,2))

ans = 0.0110

As p gets smaller, the expected value of the product becomes smaller exponentially, so it stands to reason that the number of experiments will need to increase to illustrate the result

rng(1);
p = 0.1;
experiments = 100000;
psi = rand(experiments,5) < p;
prod(mean(psi,1))

ans = 1.0246e-05

mean(prod(psi,2))

ans = 2.0000e-05

rng(1);
p = 0.1;
experiments = 100000*10; % increase experiments
psi = rand(experiments,5) < p;
prod(mean(psi,1))

ans = 1.0070e-05

mean(prod(psi,2))

ans = 1.3000e-05

[Georgios Apostolakis]

Thank you very much! I think I've made it clear on my mind now.