N=10000;
EspX=0;
VarX=0;
Esp2X=0;
EspY=0;
VarY=0;
Esp2Y=0;
S=2;
sig=2;
mu=10;

for i=1:N;
    S=2;
    while (S > 1)
        u1=rand();
        u2=rand();
        v1=2*u1-1;
        v2=2*u2-1;
        S= v1*v1 + v2*v2;
    end
    //sint=v2/Rb;
    //cost=v1/Rb;
    
    xp=v1*sqrt(-2*log(S)/S); 
    x=xp*sig+mu;
    yp=v2*sqrt(-2*log(S)/S);
    y=yp*sig+mu;
    
    EspX = x/N + EspX;
    Esp2X = (x*x)/N + Esp2X;
    
    EspY = y/N + EspY;
    Esp2Y = (y*y)/N + Esp2Y;
end;

disp(EspX);
disp(EspY);

VarX=Esp2X-EspX^2;
VarY=Esp2Y-EspY^2;

disp(VarX);
disp(VarY);
