x = [1,2,3,4,5,6];
pk = [0.1,0.2,0.1,0.3,0.1,0.2];
esp=0;
Esp2=0;
var=0;
alph=0;
n=10000;

for i = 1:n;
    somme=0;
    u=rand(i);
    for j = 1:7;
        if (u>somme) then
            somme = somme + pk(j);
            if u <=somme then 
                alph = x(j); 
            end;
        end;
    end;
    esp = esp + (alph/n);
    Esp2 = (alph^2)/n + Esp2;
end;
var = Esp2 - esp^2;
disp(var);
disp(esp);

