A=[0.4 -0.3733 ; 0.06 0.6];
B=[-0.8 -0.1867 ; 0.1371 0.8];
a=[0.3533 ; 0];
b=[1.1 ; 0.1];
N=10000;
x=grand(2,1,'def');
p=0.5;

for j=1:N
        x=bool2s(grand(1,1,'def')>p)*(A*x+a-B*x-b)+B*x+b;
        plot2d(x(1),x(2),0,rect=[0,0,1.4,0.9]);
end