clear();
X=[0 1 2 3 4]
P=[0.01 0.6 0.2 0.09 0.1]
N=10000;
Esp=0;
Var=0;
alpha=0;
Variance=0;
r=0;
cumul=0;
for i=1:N
r=rand();
J=1;
cumul=P(J);
while cumul<r
J=J+1;
cumul=cumul+ P(J);
end;
Esp=Esp + (X(J)/N);
Var=Var + (X(J) * X(J) /N);
end;
Variance = Var - (Esp*Esp);
disp(Esp);
disp(Variance);
