
// TP2 L3 proba/stats 2017-2018

// Exercice 1

//1)

p1=cdfbin("PQ",2,30,1/10,9/10)-cdfbin("PQ",1,30,1/10,9/10)
p2=cdfbin("PQ",4,30,1/10,9/10)

// Dire que Poisson(np) (pas écrit dans l'énoncé!)

// comme n>=30, p<= 0.1 et np<= 10, on peut approcher la loi binomiale par la loi de Poisson P(3) car 3=np=30*(1/10) ('TCL pas terrible')

pois1=cdfpoi("PQ",2,3)-cdfpoi("PQ",1,3)
pois2=cdfpoi("PQ",4,3)

//2) 

n=30
p=1/5
p1=cdfbin("PQ",6,n,p,1-p)-cdfbin("PQ",1,n,p,1-p)
p2=1-cdfbin("PQ",4,n,p,1-p)

// comme n>=30, np>= 5 et n(1-p)>= 5, on peut approcher la loi binomiale par la loi normale N(np,np(1-p)) ('TCL OK'))

norm1=cdfnor("PQ",6,n*p,sqrt(n*p*(1-p)))-cdfnor("PQ",1,n*p,sqrt(n*p*(1-p)))
norm2=1-cdfnor("PQ",4,n*p,sqrt(n*p*(1-p)))

// avec la correction en continuité

norm12=cdfnor("PQ",6+1/2,n*p,sqrt(n*p*(1-p)))-cdfnor("PQ",2-1/2,n*p,sqrt(n*p*(1-p)))
norm22=1-cdfnor("PQ",5-1/2,n*p,sqrt(n*p*(1-p)))

// Exercice 2

//1)

n=50
m=5
v=2
echantillon=grand(n,1,"nor",m,sqrt(v));

niveau=0.95 //0.98
alpha=1-niveau
z=cdfnor("X",0,1,1-alpha/2,alpha/2)

a=mean(echantillon)-z*stdev(echantillon)/sqrt(n)
b=mean(echantillon)+z*stdev(echantillon)/sqrt(n)
IC=[a, b]

// on veut que z*stdev(echantillon)/sqrt(n) soit inférieur à 0.1

n0=ceil( (z*stdev(echantillon)/0.01)^2 )

//2)

// Préciser à tout le monde que theta=2 par souci de cohérence.

theta=2

N=100
ech=grand(N,1,"unf",0,theta);

thetachap1=zeros(N,1);
thetachap2=zeros(N,1);

for (i = 1:N)
    thetachap1(i)=max(ech(1:i));
    thetachap2(i)=2*mean(ech(1:i));
end

plot(1:N,[thetachap1 thetachap2 ones(N,1)*theta])

// le biais de thetachap2 vaut 0 (facile)

// le biais de thetachap1 vaut theta*n/(n+1) (comme en TD: fonction de repartition 
// du max)

Nsimu=200
n=10:10:100
biais=zeros(n)

// On simule 200 échantillons de taille 10, 20...100.
// D'où la double boucle.

for i = 1:length(n)
    thetachap1=zeros(1,Nsimu);
    for j = 1:Nsimu
        ech=grand(n(i),1,"unf",0,theta);
        thetachap1(j)=max(ech);
    end
    biais(i)=mean(thetachap1)-theta
end

plot(1:100, theta*(1:100)./((1:100)+1)-theta,'green')
plot(n,biais)

// erreur quadratique moyenne

Nsimu=200
n=10:10:100
biais=zeros(n)
eqm1=zeros(n)
eqm2=zeros(n)

for i = 1:length(n)
    thetachap1=zeros(1,Nsimu);
    thetachap2=zeros(1,Nsimu);
    for j = 1:Nsimu
        ech=grand(n(i),1,"unf",0,theta);
        thetachap1(j)=max(ech);
        thetachap2(j)=2*mean(ech);
    end
    biais(i)=mean(thetachap1)-theta
    eqm1(i)=mean( (thetachap1-biais(i))^2 )
    eqm2(i)=mean( (thetachap2-theta)^2 )
end

plot(1:100, theta*(1:100)./((1:100)+1)-theta,'green')
plot(n,biais)

plot(n,[eqm1; eqm2])
legend('theta chap 1','theta chap 2');

// Exercice 3

// 1)

N=10^6
x=grand(N,1,"unf",-10,10);
y=grand(N,1,"unf",-10,10);

plot(x,y,'.')
D=(y.^3-x.*y+2*x.^2.*exp(-x.^2/100).*cos(y)<1); 
size(D)
plot(x(D),y(D),'.red')
// On sélectionne les couples qui satisfont la condition.

prop=sum(D)/N
aire=400*prop


IC=400*[prop-1.96*sqrt(prop*(1-prop)/sqrt(N)) , prop+1.96*sqrt(prop*(1-prop)/sqrt(N))]
// On modélise cela comme une v.a. de Bernouilli (dans la mesure où seulement 2 résultats possibles: soit on est dans le domaine, soit on n'est pas dedans).

// 2)

d=5
volume= %pi^(d/2)/gamma(d/2+1)

N=10^7;
x=grand(N,d,"unf",-1,1);

vol_estime=sum( sum(x.^2,2)<1 )/N*2^d 

// on peut simuler davantage de points proches de 0 car l'hypersphère est "petite". 
// Le résultat est bien plus précis.

d=10
volume= %pi^(d/2)/gamma(d/2+1)

N=10^6;
sig=0.4;
x=grand(N,d,"unf",-1,1);
y=grand(N,d,"nor",0,sig);

// On modélise les N coordonnées à la fois. On fait une loi normale avec 
// écart-type 0.4.

prop1=sum( sum(x.^2,2)<1 )/N;
prop2=mean( bool2s(sum(y.^2,2)<1) * (2*%pi*sig^2)^(d/2) .* exp(sum(y.^2,2)/2/sig^2));


vol_estime1=prop1*2^d 
vol_estime2=prop2 


//Notons 1_s l'indicatrice de la sphère. Notons f la densité de la loi N(0,sigma).
//
//Alors E[1_s(U_1,...,U_d)] = \int_s  1/2^d du_1...du_d
//
//et grâce à la loi des grands nombres, on estime l'espérance par la moyenne empirique de N réalisations de (U_1,...,U_d). D'où la formule : volume = 2^d * moyenne empirique des 1_s(U_1,...,U_d)
//
//On peut aussi écrire :
//
//volume = \int_[-1,1]^d  1_s(u_1,...,u_d)  du_1...du_d
//
//            = \int_R^d  1_s(u_1,...,u_d)/f(u_1,...,u_d) * f(u_1,...,u_d) du_1...du_d
//
//            = E[ 1_s(Z_1,...,Z_d)/f(Z_1,...,Z_d) ] où les Z_i sont de loi normale N(0,sigma).
// (en effet, f(u_1,...,u_d) du_1...du_d est une densité de probabilité).
//
//et idem, on estime l'espérance par la moyenne empirique. Cela s'appelle de l'échantillonnage préférentiel (on échantillonne là où il se passe quelque chose).

// Exercice 5

// une partie

c=20
a=10
A=a

while (a>0)&(a<c)
    jeu=2*grand(1,1,"bin",1,0.5)-1;
    a=a+jeu;
    A=[A, a];
end

plot2d2(1:length(A),A)

// quelle est la proba de gagner avec a=10 ?

Nsimu=10^3
c=20

vainqueur=zeros(Nsimu,1);

for i=1:Nsimu
    a=10;
    while (a>0)&(a<c)
        a=a+2*grand(1,1,"bin",1,0.5)-1;
    end
    if (a==c) 
        vainqueur(i)=1;
    end
end

proba=mean(vainqueur)

// pour tout a

Nsimu=10^3
c=20
fortune_init=0:c;

proba=zeros(fortune_init);

for j=1:length(fortune_init)
    vainqueur=zeros(Nsimu,1);
    for i=1:Nsimu
        a=fortune_init(j);
        while (a>0)&(a<c)
            a=a+2*grand(1,1,"bin",1,0.5)-1;
        end
        if (a==c) 
            vainqueur(i)=1;
        end
        proba(j)=mean(vainqueur);
    end
end

plot(fortune_init,proba)



