% tp1 M1 info DPAM
%
% sol. FC

clear all
close all
clc

set(0,'defaultaxesfontsize',15); %option



%% 1. Loi Binomiale de paramètre p :  Bin(x;n,p) %%

% cas 1 : n=30, p=0.5

n1=30; 
p1 = 0.5;
x1 = 0:50;
f1 = binopdf(x1,n1,p1);

figure;
plot(x1,f1,'r.');
% legend(['Bin(',int2str(n1),',', num2str(p1),')']);

% cas 2 : n=30, p=0.7

n2=30; 
p2=0.7; 
x2 = 0:50;
f2 = binopdf(x2,n2,p2);

hold on; plot(x2,f2,'g.');
% legend(['Bin(',int2str(n2),',', num2str(p2),')']);
% cas 3 : n=50, p=0.4
n3 = 50; p3 = 0.4;
x3 = 0:50;
f3 = binopdf(x3,n3,p3);

hold on; plot(x3,f3,'b.');
legend(['Bin(',int2str(n1),',', num2str(p1),')'], ['Bin(',int2str(n2),',', num2str(p2),')'], ['Bin(',int2str(n3),',', num2str(p3),')']);
xlabel('x'); ylabel('p(x)')
title('Loi Binomiale : Bin(n,p)')


%% Loi Normale Univariée  N(x;mu,sigma)%%

% 1. loi centrée réduite : mu = 0, sigma = 1
x = -4:0.01:7;

mu1 = 0; sigma1 = 1;
f1 = normpdf(x,mu1,sigma1); 


figure,
plot(x,f1,'r')

% 2. loi normale : mu = 0.5, sigma = 1
mu2 = 2; sigma2 = 1.5;
f2 = normpdf(x,mu2,sigma2);

hold on; plot(x,f2,'g')


% 3. loi normale : mu = 0.5, sigma = 0.3
mu3 = mu2; sigma3 = 0.6 ;
f3 = normpdf(x,mu3,sigma3);

hold on; plot(x,f3,'b')
legend(['N(',int2str(mu1),',', num2str(sigma1),'^2)'], ...
    ['N(',int2str(mu2),',', num2str(sigma2),'^2)'],...
    ['N(',int2str(mu3),',', num2str(sigma3),'^2)']);
xlabel('x'); ylabel('f(x)')
title('Loi Normale : N(\mu,\sigma^2)')
xlim([-4 7]);

%% Application de la loi Normale multivariée N(x;Mu,Sigma)%%(cas 2d)

mu = zeros(1,2);
sigma = eye(2);

% axis
x1liminf= -5; x1limsup = 5;
x2liminf= -5; x2limsup = 5;
nx = 100; ny = 100;
ngrid=100;
%[X1,X2] = meshgrid(linspace(x1liminf,x1limsup,nx),linspace(x2liminf,x2limsup,ny));
[X1,X2] = meshgrid(linspace(-5,5,100),linspace(-5,5,100));

X1 = X1(:); X2 = X2(:);
X  = [X1 X2];

f = mvnpdf(X,mu,sigma);

%% commandes : contours, mesh et surf pour afficher le fdp
X1=reshape(X1,ngrid,ngrid);
X2=reshape(X2,ngrid,ngrid);
f=reshape(f,ngrid,ngrid);

%% contour
figure, contour(X1,X2,f); 
xlabel('x1');ylabel('x2');title(['contour de f(x) = N(0, I)']);

%
niveau = 0.01;
figure, contour(X1,X2,f,niveau);
xlabel('x1');ylabel('x2');title(['contour de f(x) à un niveau ',num2str(niveau)]);


%% mesh
C=.8*ones(ngrid,ngrid);

figure, mesh(X1,X2,f,C);
xlabel('x1');ylabel('x2');zlabel('f(x)')

%% surf
figure, surf(X1,X2,f);
xlabel('x1');ylabel('x2');zlabel('f(x)')

%scatter(X(:,1),X(:,2),'r.');


%% Simulation d'une Loi Normale Multivariée N(x;Mu,Sigma)%%
% 1.  Loi Normale Multivariée centrée réduite
n = 500; % nombre de point

mu1 = zeros(1,2);
sigma1 = eye(2);

Y = mvnrnd(mu1, sigma1,n);

f1 = mvnpdf(X,mu1,sigma1);

figure;
scatter(Y(:,1),Y(:,2),’r.’);
hold on;
plot(mu1(1),mu1(2),'k.', 'linewidth',2,'markersize',20);
hold on;
draw_ellipse(mu1,sigma1, 'b');


% 2. 
mu2 = [-4 1.5]; 
sigma2 = [0.3 0.1;0.1 0.15]; 

Y = mvnrnd(mu2, sigma2,n);
f2 = mvnpdf(X,mu2,sigma2);

hold on
scatter(Y(:,1),Y(:,2),’g.’);
hold on;
plot(mu2(1),mu2(2),'k.', 'linewidth',2,'markersize',20);
hold on;
draw_ellipse(mu2,sigma2, 'b');


% 3. 
mu3 = [5 1]; 
sigma3 = [1 -.3;-.3 0.2];

Y = mvnrnd(mu3, sigma3,n);
f3 = mvnpdf(X,mu3,sigma3);

hold on
scatter(Y(:,1),Y(:,2),’m.’);
hold on;
plot(mu3(1),mu3(2),'k.', 'linewidth',2,'markersize',20);
hold on;
draw_ellipse(mu3,sigma3, 'b');

xlabel('x1');ylabel('x2')
xlim([-6 8])
box on


%% loi normale univariée : histogramme et densité

n = 200;
mu = 0; 
sigma = 1;

x = linspace(-4,4,100);% -4:0.1:4;
%y = randn(n,1);
y = normrnd(mu,sigma,n,1); 

f1 = normpdf(x,mu,sigma); 

% figure,
% hist(y,x,15);
% hold on; plot(x,f1,'r','linewidth',2)
% xlim([-4 4]);
% pause

figure;
normdenhist(y,mu, sigma);
xlim([-4 4]);

%hold on, plot(x,zeros(size(x)),'kx')


%% Estimation de densité : cas de la loi normale (univariée)

% vrais paramètres
mu = 0; sigma = 1; 


% taille respectives des échantillons
n1 = 20; n2 = 80; n3 = 300;
x1 = normrnd(mu,sigma,n1,1);
x2 = normrnd(mu,sigma,n2,1);
x3 = normrnd(mu,sigma,n3,1);

% Estimation de mu et sigma^2
m1 = mean(x1); s1 = var(x1);
m2 = mean(x2); s2 = var(x2);
m3 = mean(x3); s3 = var(x3);

s1 = sqrt(s1);
s2 = sqrt(s2);
s3 = sqrt(s3);

% ou 
% s1 = sqrt((1/(n1-1))*(x1 - repmat(m1,n1,1))'*(x1 - repmat(m1,n1,1)));
% s2 = sqrt((1/(n2-1))*(x2 - repmat(m2,n2,1))'*(x2 - repmat(m2,n2,1)));
% s3 = sqrt((1/(n3-1))*(x3 - repmat(m3,n3,1))'*(x3 - repmat(m3,n3,1)));
 
x = -4:0.1:4;
% vrai densité
f = normpdf(x,mu,sigma); 

% densités estimées
f1 = normpdf(x,m1,s1); 
f2 = normpdf(x,m2,s2); 
f3 = normpdf(x,m3,s3); 

%plots
figure, plot(x,f,'b','linewidth',2)
legend('vraie densité')
hold on; plot(x,f1,'r')
hold on; plot(x,f2,'g')
hold on; plot(x,f3,'m')
legend('vraie densité', ...
    ['n = ', num2str(n1)],...
    ['n = ', num2str(n2)],...
    ['n = ', num2str(n3)]);
xlabel('x'); ylabel('f(x)');

%% cdf

% vrai cdf
cdf = normcdf(x,mu,sigma); 

% densités estimées
cdf1 = normcdf(x,m1,s1); 
cdf2 = normcdf(x,m2,s2); 
cdf3 = normcdf(x,m3,s3); 

%plots
figure, plot(x,cdf,'b','linewidth',2)
legend('vraie densité')
hold on; plot(x,cdf1,'r')
hold on; plot(x,cdf2,'g')
hold on; plot(x,cdf3,'m')
legend('vraie cdf', ...
    ['n = ', num2str(n1)],...
    ['n = ', num2str(n2)],...
    ['n = ', num2str(n3)]);
xlabel('x'); ylabel('cdf(x)');


%% Estimation de densité : cas de la loi exponentielle

% vrais paramètre
lambda = 1.5;


% taille respectives des échantillons
n1 = 20; n2 = 50; n3 = 100;
x1 = exprnd(lambda,n1,1);
x2 = exprnd(lambda,n2,1);
x3 = exprnd(lambda,n3,1);

% Estimation de lambda (MV)
lambda1 = mean(x1);  
lambda2 = mean(x2);  
lambda3 = mean(x3); 


x = 0:0.1:8;
% vrai densité
f = exppdf(x,lambda); 

% densités estimées
f1 = exppdf(x,lambda1); 
f2 = exppdf(x,lambda2); 
f3 = exppdf(x,lambda3); 

%plots
figure, plot(x,f,'b','linewidth',2)
legend('vraie densité')
hold on; plot(x,f1,'r')
hold on; plot(x,f2,'g')
hold on; plot(x,f3,'m')
legend('vraie densité', ...
    ['n = ', num2str(n1)],...
    ['n = ', num2str(n2)],...
    ['n = ', num2str(n3)]);
xlabel('x'); ylabel('f(x)');
title(['loi exponentielle : \lambda = ',num2str(lambda)]);
%% cdf

% vrai cdf
cdf = expcdf(x,lambda); 

% cdf estimées
cdf1 = expcdf(x,lambda1); 
cdf2 = expcdf(x,lambda2); 
cdf3 = expcdf(x,lambda3); 

%plots
figure, plot(x,cdf,'b','linewidth',2)
legend('vraie densité')
hold on; plot(x,cdf1,'r')
hold on; plot(x,cdf2,'g')
hold on; plot(x,cdf3,'m')
legend('vraie cdf', ...
    ['n = ', num2str(n1)],...
    ['n = ', num2str(n2)],...
    ['n = ', num2str(n3)]);
xlabel('x'); ylabel('cdf(x)');
