function y=couleursroulette(n)
   x=grand(1,n,"uin",1,38); 
   y=1*[x<=18]+2*[x>20];  
endfunction

function z=gainsroulette(n)
    x=grand(1,n,"uin",1,38);
    y=1*[x<=18]-1*[x>18]; 
    y=[0 y]; 
    z=cumsum(y); 
    plot(0:n,z); 
endfunction

function p=estimeprobaroulette(k,n)
    p=0;
    for i=1:n do 
        x=grand(1,k,"uin",1,38);
    y=1*[x<=18]-1*[x>18]; 
    p=p+[sum(y)>=0]/n; 
end 
endfunction

function f=estimeprobaE1(n); 
    f=0;
    for i=1:n do 
        x=rand(1:10^7); 
        f=f+sum([x <= 0.5])/(n*10^7); 
    end
endfunction





function y=f1(x) 
    y=x^10
endfunction

function y=f2(x) 
    if x==0 then y=0; 
        else y=1/x; 
    end
endfunction


function e=estimeroulette(n) 
    x=rand(n,3); 
    y=[x<1/19]; // on teste si le lancer 
    z=[x>=10/19]; 
    r=[sum(y,'c')==1]; 
    s=[sum(z,'c')==1]; 
    t=[r+s==2];
    e=sum(t)/n; 
endfunction

function r=rejet1(n) 
 x = grand(1:n,"poi",22); 
y=and([x>=25; x<=40],'r');
z=find(y);  
r=x(z);
endfunction

function z=rejet2(k) 
      c=1; z=zeros(1:k);  
      while c <= k do 
          x=grand(1,1,"poi",22); 
          if and([x >=25 x<= 40]) then z(c)=x; c=c+1; 
          end
      end  
endfunction

function s=rejet3(k)
    l=0;
    while l<k do 
        n=(1.1)*k/0.28; 
        x=grand(1:n,"poi",22); 
        y=and([x>=25; x<=40],'r');
        z=find(y);
        l=length(z);  
    end
   r=x(z);
   s=r(1:k); 
endfunction

function s=rejet4(k)
    r=rejet1(1.1*k/0.28);
    while length(r)<k do r=rejet1(1.1*k/0.28); end 
    s=r(1:k);
endfunction

function b=ber(n,m,p)
   x= rand(n,m);
   b=1.*[x<p] 
endfunction

function s=bin(m,n,p)
s=sum(ber(n,m,p),'r')
endfunction

function u=unifo(m,n)
    x=rand(1:m);
    u=ceil(n.*x);
endfunction

function g=geo(p) 
    g=1; 
    while rand()>p do g=g+1; 
    end 
endfunction

function g=geom(n,p) 
    g=ones(1:n); 
    for i=1:n do 
     while rand()>p do g(i)=g(i)+1; 
     end 
    end
    g
endfunction

function g=geo2(p)
     x=rand();
     a=p; g=1; 
    while x>a do  g=g+1;   a=1-(1-p)^g;
    end
    g
endfunction


function g=geom2(n,p)
   x=rand(1:n); 
   m=max(x); a=p; g=ones(1:n); k=1; 
    while m>a do g=g+[x>a]; k=k+1; a=1-(1-p)^k; end 
    g 
endfunction

function k=poiss(lambda) 
    x=rand(); a=%e^(-lambda); k=0; 
    while x>a do k=k+1; a=a+%e^(-lambda)*lambda^k/factorial(k); 
    end
    k
endfunction

function y=poisson(n,lambda)
   x=rand(1:n); m=max(x); 
   a=%e^(-lambda); y=zeros(1:n); k=0; 
  while m>a do 
      y=y+[x>a]; 
      k=k+1; 
      a=a+%e^(-lambda)*(lambda^k)/factorial(k);  
      end 
  y
endfunction

function h=poursimulerunif(u)
    h=%pi+2*%pi*u; 
endfunction

function h=poursimulerexp1(u)
    h=-log(1-u);
endfunction

function p=pourverifdensiteexp(x)
    p=%e^(-x);
endfunction

function h=poursimulexp(lambda,u)
      h= -(1/lambda)*log(1-u); 
endfunction

function h=poursimularcsine(u) 
    h=1/2-cos(%pi*u)/2; 
endfunction

function c=comparedensityarcsine(x)
    if x <=0 then c=0; 
        else if x>=1 then c=0;
    else 
    c=1/(%pi*(x*(1-x))^(1/2));
end 
end 
endfunction

function x=poursimulnor(u) 
    x=[sqrt(-2*log(u(1)))*cos(2*%pi*u(2));sqrt(-2*log(u(1)))*sin(2*%pi*u(2))]; 
endfunction


function e=expo(x) 
    e=exp(x); 
    endfunction 


function y=expon(n,lambda)
    x=rand(1:n); 
    y=feval(lambda,x,poursimulexp);
endfunction


function f=partiel(x)
    if x <= 0 then f=0; 
    else if x < 1 then f=19*x/30;         
    else c=1; f=19/30; 
        while x >=c do f=f+1/3^(c+1); 
            if c==2 then f=f+1/5; end 
            c=c+1;          
          end
      end
  end
endfunction


function histo3dnor(n) 
  M=grand(n,2,"nor",0,1);
x=[-3:0.25:3];
y=[-3:0.25:3];
for i=1:24 do
for j=1:24 do
B=[x(i) <= M(:,1) M(:,1) <= x(i+1) y(j)<= M(:,2) M(:,2) <=y(j+1) ];
C=and(B,'c');
D(i,j) = sum(C);
end
end
hist3d(D, flag=[2 2 0]);
endfunction


function histo3dexp(n) 
  M=grand(n,2,"exp",1);
x=[-3:0.25:3];
y=[-3:0.25:3];
for i=1:24 do
for j=1:24 do
B=[x(i) <= M(:,1) M(:,1) <= x(i+1) y(j)<= M(:,2) M(:,2) <=y(j+1) ];
C=and(B,'c');
D(i,j) = sum(C);
end
end
hist3d(D, flag=[2 2 0]);
endfunction

function y=f2(x) 
    if and([x>0 x<1]) then y=1/(2*x^(1/2));
        else y=0; 
    end
endfunction

function i=veriff2density 
    i = abs(intg(0,1,f2)-1)<10^(-5)
endfunction

function y=repartf2(x)
    if x<0 then y=0; 
        else y=intg(0,x,f2); 
    end
endfunction

function y=reparttheorique(x) 
    if x<0 then y=0; 
    else if x<1 then y=x^(1/2); 
        else y=1;
            end 
    end
endfunction

function y=forcomparerepartf2
    a=[0:0.01:1]; 
    for i=1:101 do 
        y(i)=intg(0,a(i),f2);
    end
endfunction 

function e=qqsesp(n) 
     u=grand(n,1,"unf",0,1);
    x=u.^2;
   e(1)=sum(x)/n; e(2)=sum(x.^2)/n; e(3)=e(2)-e(1)^2;   
endfunction

function representplane(n)
    x=grand(n,2,"unf",0,1);
    y=x.^2; 
    for i=1:n do 
        plot(y(i,1),y(i,2),'x');
    end
endfunction


function hist3dcouple(n)
     x=grand(n,2,"unf",0,1); 
     y=x.^2;  
a=[-0.1:0.1:1.1]; 
b=[-0.1:0.1:1.1]; 
for i=1:12 do 
    for j=1:12 do 
        B=[a(i) <=y(:,1) y(:,1)<=a(i+1) b(j)<=y(:,2) y(:,2) <= b(j+1)]; 
        C=and(B,'c'); 
        D(i,j)=sum(C);
    end
end
hist3d(D);
endfunction 


function z=conditionnement1(n) 
    x=grand(n,2,"unf",0,1); 
     y=x.^2;  
     c=1; 
     for i=1:n do 
         if y(i,1)<y(i,2) then z(c,1)=y(i,1); z(c,2)=y(i,2); c=c+1; end 
     end
endfunction


function hist3dcouple2(n)
     y=conditionnement1(n);
a=[-0.1:0.1:1.1]; 
b=[-0.1:0.1:1.1]; 
for i=1:12 do 
    for j=1:12 do 
        B=[a(i) <=y(:,1) y(:,1)<=a(i+1) b(j)<=y(:,2) y(:,2) <= b(j+1)]; 
        C=and(B,'c'); 
        D(i,j)=sum(C);
    end
end
hist3d(D);
endfunction 


function hist3dnor(n)
    M= grand(n,2,"nor",0,1); 
x=[-3:0.25:3]; 
y=[-3:0.25:3]; 
for i=1:24 do 
    for j=1:24 do 
        B=[x(i) <=M(:,1) M(:,1)<=x(i+1) y(j)<=M(:,2) M(:,2) <= y(j+1)]; 
        C=and(B,'c'); 
        D(i,j)=sum(C)/sqrt(n);
    end
end
hist3d(D);
endfunction 

function norber(n)
    x=grand(n,1,"nor",0,1); 
    eps=grand(n,1,"uin",0,1); 
    z=2 .* eps -1; 
    y=z .* x; 
    histplot([-4:0.1:4],y); 
    u=[-3:0.25:3]; 
    v=[-3:0.25:3]; 
for i=1:24 do 
    for j=1:24 do 
        B=[u(i)<=x u(i+1)>x  v(j)<=y  y < v(j+1)]; 
        C=and(B,'c'); 
        D(i,j)=sum(C)/sqrt(n);
    end
end
// clf()
hist3d(D);
endfunction

function z=unifdisque(n) 
    x=grand(2,n,"unf",-1,1); 
    c=1; 
    for i=1:n do
    if x(1,i)^2+x(2,i)^2 <= 1 then z(c,1)=x(1,i); z(c,2)=x(2,i); c=c+1; end 
    end
endfunction

function z=betterunifdisque(n) 
    x=grand(2,n,"unf",-1,1); 
    y=x.^2; s=sum(y,'r'); 
    f=find(s <=1); 
    z=x(:,f)';
endfunction

function hist3dunifdisque(n)
    z=betterunifdisque(n); 
    x=[-1:0.1:1]; 
y=[-1:0.1:1]; 
for i=1:20 do 
    for j=1:20 do 
        B=[x(i) <=z(:,1) z(:,1)<=x(i+1) y(j)<=z(:,2) z(:,2) <= y(j+1)]; 
        C=and(B,'c'); 
        D(i,j)=sum(C)/sqrt(n);
    end
end
hist3d(D);
endfunction


function x=simulgaussienne(m,S,n);
    d=length(m);
         y=grand(d,n,"nor",0,1); 
     [P,D]=schur(S); 
     x=P'*sqrt(D)*y+m*ones(1,n); 
endfunction

function hist3dgauss(m,S,n)
    M=simulgaussienne(m,S,n)'; 
x=[-3:0.25:3]; 
y=[-3:0.25:3]; 
for i=1:24 do 
    for j=1:24 do 
        B=[x(i) <=M(:,1) M(:,1)<x(i+1) y(j)<=M(:,2) M(:,2) < y(j+1)]; 
        C=and(B,'c'); 
        D(i,j)=sum(C)/sqrt(n);
    end
end
hist3d(D);
endfunction 

function repartgeo(p)
    x=grand(1,10^6,"geom",p); 
    a=[0:0.05/p:5/p];
     n=length(a); 
     y=zeros(1:n); 
    for i=1:n do 
        y(i)=sum([x <= a(i)])/10^6; 
    end
    plot(a,y); 
endfunction

function z=repartexp(x)
    z=1-%e^(-x); 
endfunction

function repartgeobis(p)
    x=p*grand(1,10^6,"geom",p); 
    a=[0:0.05:5];
     n=length(a); 
     y=zeros(1:n); 
    for i=1:n do 
        y(i)=sum([x <= a(i)])/10^6; 
    end
    z=feval(a,repartexp);
    plot(a,y); 
    plot2d(a,z,style=-2);
endfunction

function histosgeomexp(n,p)
    x=grand(1,n,"exp",1);
    y=p*grand(1,n,"geom",p); 
    a=[0:0.1:10];
    histplot(a,x);
    histplot(a,y,style=2);
endfunction

function s=forLGNunf(n,m,a,b)
    s=[]; 
    for i=1:n do 
        x=grand(m,1000,"unf",a,b); 
        y=sum(x,'r')/m; 
        s=[s y];
    end
    disp([min(s) max(s)]); 
    a=[(a+b)/2-1:0.001:(a+b)/2+1];
    histplot(a,s); 
endfunction

function s=forLGNexp(n,m,lambda)
    s=[]; 
    for i=1:n do 
        x=grand(m,1000,"exp",lambda); 
        y=sum(x,'r')/m; 
        s=[s y];
    end
    disp([min(s) max(s)]); 
    a=[lambda-1:0.001:lambda+1];
    histplot(a,s); 
endfunction

function s=forLGNnor(n,m,mu,v)
    s=[]; 
    for i=1:n do 
        x=grand(m,1000,"nor",mu,v); 
        y=sum(x,'r')/m; 
        s=[s y];
    end
    disp([min(s) max(s)]); 
    a=[mu-1:0.001:mu+1];
    histplot(a,s); 
endfunction

function s=forLGNpoi(n,m,lambda)
    s=[]; 
    for i=1:n do 
        x=grand(m,1000,"nor",lambda); 
        y=sum(x,'r')/m; 
        s=[s y];
    end
    disp([min(s) max(s)]); 
    a=[lambda-1:0.001:lambda+1];
    histplot(a,s); 
endfunction

function y=forcauchy(x)
    y=tan(%pi*(x-1/2));
endfunction

function cexcauchy(n,m)
    s=[]; 
    for i=1:n do 
        u=grand(m,1000,"unf",0,1);
        x=feval(u,forcauchy);  
        y=sum(x,'r')/m; 
        s=[s y];
    end
    disp([min(s) max(s)]); 
    a=[-10:0.5:10];
    histplot(a,s);
endfunction

function y=cauchydensity(x)
    y=1/(%pi*(1+x^2)); 
endfunction

function g=gauss(x)
    g=%e^(-x^2/2)/(sqrt(2*%pi)); 
endfunction

function TCL1(n,m)
    stacksize('max');  
    s=zeros(1:n)
    for i=1:n do 
        x=2*grand(1,m,"uin",0,1)-1; 
        s(i)=sum(x)/sqrt(m); 
    end
    histplot(100,s);  
    a=[-4:0.01:4]; 
    y=feval(a,gauss);
    plot(a,y);
endfunction

function s=forTCL(n,m)
   s=[];
    for i=1:n do 
        x=2*grand(m,1000,"uin",0,1)-1; 
        y=sum(x,'r')/sqrt(m); 
        s=[s y];
    end
endfunction

function s=forTCL2(n,m)
    s=[];
    for i=1:n do 
        x=grand(m,1000,"exp",1)-1; 
        y=sum(x,'r')/sqrt(m); 
        s=[s y];
    end
endfunction

    
function s=forTCL3(n,m)
    s=[];
    for i=1:n do 
        x=grand(m,1000,"gam",3,1)-3; 
        y=sum(x,'r')/sqrt(3*m); 
        s=[s y];
    end
endfunction    
    


function TCL(s)
    histplot(100,s);  
    a=[-4:0.01:4]; 
    y=feval(a,gauss);
    plot(a,y);
endfunction
