% Lezione di lab TTR  -9/10/2013- Legge di Bayes - Marco Cristani
close all;
clear all;
clc
% ---------------
% Esercizio 1
% ---------------
fprintf(' ---------------\n')
fprintf(' Esercizio 1\n')
fprintf(' ---------------\n')
%-- Visualizzazione campioni
mu1         =   1;
sigma1      =   2;

mu2         =   20;
sigma2      =   3;

num_s       =   500;

%fprintf('Generazione campioni Gaussiani 1D...');
[Sample_set1] =   normrnd(mu1,sigma1,num_s,1);
[Sample_set2] =   normrnd(mu2,sigma2,num_s,1);
%fprintf('OK.\n');

%Come fare intuitivamente a capire che sono distribuiti in modo normale?
%Con l'istogramma
[histo1,c1]=hist( Sample_set1,32);
[histo2,c2]=hist( Sample_set2,32);

h0=figure;
set(gcf,'name','visualizzazione istogrammi')
subplot(211); bar(c1,histo1);
subplot(212); bar(c2,histo2);


Data        =   [Sample_set1;Sample_set2];

%-- Plot campioni

maxData     =   max(Data);
minData     =   min(Data);


h_G1d = figure; hold on;
h_S1 = scatter(Sample_set1,zeros(num_s,1),10,'b');
h_S2 = scatter(Sample_set2,zeros(num_s,1),10,'r');



axis([minData-1,maxData+1,0,2])
set(gcf,'name','Campioni Gaussiani')
set(gca,'ytick',[]);
h_m1 = line([mu1;mu1],[0;2]); 
set(h_m1,'color',[0,0,1]);
text(mu1+0.05,1,'\mu1','color',[0,0,1])

h_m2 = line([mu2;mu2],[0;2]);
set(h_m2,'color',[1,0,0]);
text(mu2+0.05,1,'\mu2','color',[1,0,0])

h_leg = legend([h_S1(1);h_S2(1);h_m1;h_m2],{strcat('N1: \mu1=',num2str(mu1),' \sigma1=',num2str(sigma1));...
    ,strcat('N2: \mu2=',num2str(mu2),' \sigma2=',num2str(sigma2));...
    '\mu1';...
    '\mu2'});


fprintf('mu1 = %3.5f; dev-std1 =%3.5f\n',mu1,sigma1)
fprintf('mu1 = %3.5f; dev-std1 =%3.5f\n',mu2,sigma2)

%-- Stima media e varianza

%fprintf('Stima media e varianza...')
est_mu1 = mean(Sample_set1);
est_sigma1 = std(Sample_set1);

est_mu2 = mean(Sample_set2);
est_sigma2 = std(Sample_set2);
%fprintf('OK.\n');


fprintf('mu1 stimata = %3.5f; dev-std1 stimata =%3.5f\n',est_mu1,est_sigma1)
fprintf('mu1 stimata = %3.5f; dev-std1 stimata =%3.5f\n',est_mu2,est_sigma2)

%-- Plot stime
set(h_G1d,'name','Campioni Gaussiani, confronto medie reali e stimate')
h_estm1 = line([est_mu1;est_mu1],[0;2]); 
set(h_estm1,'color',[0,0,1],'LineStyle','--');
text(est_mu1+0.05,1.5,'est. \mu1','color',[0,0,1])

h_estm2 = line([est_mu2;est_mu2],[0;2]);
set(h_estm2,'color',[1,0,0],'LineStyle','--');
text(est_mu2+0.05,1.5,'est. \mu2','color',[1,0,0])

delete(h_leg)

h_leg = legend([h_S1(1);h_S2(1);h_m1;h_m2;h_estm1;h_estm2],{strcat('N1: \mu1=',num2str(mu1),' \sigma1=',num2str(sigma1));...
    ,strcat('N2: \mu2=',num2str(mu2),' \sigma2=',num2str(sigma2));...
    ' \mu1';...
    ' \mu2';...
    'est. \mu1';...
    'est. \mu2'});

%-- Errore stime
%fprintf('Calcolo errore stime...');
err_mu1 = abs(est_mu1-mu1);
err_sigma1 = abs(est_sigma1-sigma1);

err_mu2 = abs(est_mu2-mu2);
err_sigma2 = abs(est_sigma2-sigma2);
%fprintf('OK \n');

fprintf('Errore mu1 = %3.5f; Errore dev-std1 =%3.5f\n',err_mu1,err_sigma1);
fprintf('Errore mu2 = %3.5f; Errore dev-std2 =%3.5f\n',err_mu2,err_sigma2);

fprintf('FINE! Premi un tasto per continuare!\n')
pause

% ---------------
% Esercizio 2
% ---------------
fprintf(' ---------------\n')
fprintf(' Esercizio 2\n')
fprintf(' ---------------\n')

% Data        =   [Sample_set1;Sample_set2];

% valuto la likelihood dei miei punti rispetto ad una data
% parametrizzazione Gaussiana
mu3 = 10.5;
dev_std3 = 2;

lik_ind3 = normpdf(Data,mu3,dev_std3);


h3 = figure;
scatter(Data,lik_ind3);

% la figura mi dice che mu3 e sigma3 c'azzeccano poco con i dati. Voglio
% quantificare la mia assunzione. Calcolo la likelihood.

lik_tot3 = prod(lik_ind3);
fprintf('Proba che theta_3 abbia generato i dati: %f\n',lik_tot3);

llik_tot3 = sum(log(lik_ind3));

% la likelihood risulta molto bassa... sbaglio qualcosa! provo quindi a
% valutare la likelihood che i punti siano effettivamente distribuiti
% secondo le distribuzioni che li hanno generati;


lik_ind1 = normpdf(Sample_set1,mu1,sigma1);

lik_ind2 = normpdf(Sample_set2,mu2,sigma2);

lik_tot12=[lik_ind1;lik_ind2];

llik_tot12 = sum(log(lik_tot12));

llik_tot3
llik_tot12



%Manipolazione di distribuzioni gaussiane, unimodali e multivariate
clear all;
A= normrnd(0,2,10000,2);

h1 = figure;
scatter(A(:,1),A(:,2));
axis equal

A2d = mvnrnd([0,0],[4,0;0,4],10000);
h2 = figure;
scatter(A2d(:,1),A2d(:,2));
axis equal

% Inserisco una leggera covarianza tra le variabili
A2dc = mvnrnd([0,0],[4,0.4;0.4,4],10000);
h2 = figure;
scatter(A2dc(:,1),A2dc(:,2));
axis equal
fprintf('FINE! Premi un tasto per continuare!\n')
pause

% ---------------------------------------
% Esercizio 3 - Riconoscimento di facce
% ---------------------------------------
fprintf(' ---------------\n')
fprintf(' Esercizio 3\n')
fprintf(' ---------------\n')
A=rgb2gray(imread('1.jpg'));
A=A(1:4:end,1:4:end);
rows = size(A,1);
cols = size(A,2); 

figure; 
set(gcf,'name','Medie e varianze, metodo veloce')
colormap gray; subplot 221 ; imagesc(A); 
title('Mean class 1');
Mu1 = A(:);
Sigma1 = ones(rows*cols,1).*10;
subplot 222 ; imagesc(reshape(Sigma1,rows,cols)); colorbar; 
title('Sigma class 1');

B=rgb2gray(imread('12.jpg'));
B=B(1:4:end,1:4:end);
subplot 223; imagesc(B); 
title('Mean class 2');
Mu2 = B(:);
Sigma2 = ones(rows*cols,1).*10;
subplot 224 ; imagesc(reshape(Sigma2,rows,cols)); colorbar; 
title('Sigma class 2');

C1=[];
C2=[];
for i=1:12
    test = rgb2gray(imread(strcat(num2str(i),'.jpg')));
    test = test(1:4:end,1:4:end);
    test=test(:);
    LK1 = sum(log(normpdf(double(test),double(Mu1),double(Sigma1))));
    LK2 = sum(log(normpdf(double(test),double(Mu2),double(Sigma2))));
    if LK1 >LK2
        C1 = [C1,i];
    else
        C2 = [C2,i];
    end
end
latoC1=ceil(sqrt(length(C1)));
latoC2=ceil(sqrt(length(C2)));
figure;
set(gcf,'name','Classe 1')
count=1;
for i=1:length(C1)
    
    subplot(latoC1,latoC1,i);
    test = rgb2gray(imread(strcat(num2str(C1(i)),'.jpg')));
    imagesc(test);axis off; colormap gray;
    
end

figure;
set(gcf,'name','Classe 2')
count=1;
for i=1:length(C2)
    
    subplot(latoC2,latoC2,i);
    test = rgb2gray(imread(strcat(num2str(C2(i)),'.jpg')));
    imagesc(test);axis off; colormap gray;
    
end    
fprintf('FINE! Premi un tasto per continuare!\n')
pause

