% Quarto  laboratorio 3 dicembre 2002


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Argomenti: Simulazione di variabili aleatorie, 
%            convergenza di variabili aleatorie
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1. SIMULAZIONE DI VARIABILI ALEATORIE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Generazione della uniforme continua su (0,1): 

% RAND: Uniformly distributed random numbers.  Per generare una matrice
% N*M di variabili pseudo aleatorie indipendenti e uniformi in
% (0,1)usiamo

U = rand(N,M)
u = rand(N,1) % fornisce un vettore colonna di lunghezza N.

% Una prima indagine (grafica) per testare la bontà del generatore è la
% seguente:

% primo passo: Fissiamo n ``grande'', per esempio provare con
% n = 500, 5000, 50000 e generiamo 
n=500
u = rand(n,1);

% secondo passo: Costruiamo due vettori di lunghezza 2k il primo -diciamo
% u1- costituito dagli elementi di u di posto pari e il secondo -diciamo
% u2- da quelli di posto dispari:

u1=u(1:2:n);
u2=u(2:2:n);

% A questo punto plottiamo u1 versus u2:

plot(u1,u2,'o') 

% Se i punti sono uniformemente distribuiti sul quadrato [0,1]*[0,1]
% allora possiamo congetturare che non ci sia correlazione fra i vari
% elementi del vettore u.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% Generazione di altre variabili diverse dalla uniforme. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Matlab le fornisce di default, tutte partono dal generare numeri
% indipendenti uniformi (0,1) poi trasformati.
% binornd, unifrnd, normrnd, poissrnd...

% esempio:
u = normrnd(0,1,n,1);
u1=u(1:2:n);
u2=u(2:2:n);
plot(u1,u2,'o')

% Esercizio 1. Sia U una variabile aleatoria uniforme in (0,1).  sia X
% definita come: se U<p allora X=1, se U>=p allora X=0.  Qual è la legge
% di X?
% Soluzione. X è una variabile aleatoria bernoulliana di parametro p.
% Infatti: {X=1} = {U<p} cosicchè P(X=1) = P(U<p) = p.

% Esercizio 2. Si generino direttamente N (=45) numeri dalla
% bernoulliana di parametro p (=0.3)
% Soluzione.  Il comando successivo genera un vettore colonna di 0,1,
% secondo la seguente regola: x(i) = 1 se u(i) < p.
 
rand(n,1) < p

% Esercizio 3. Si generino direttamente n (=100) numeri dalla binomiale
% di parametri (N=45,p=0.3):

% Risp   Poichè una Bin(N,p) è una somma di N bernoulli indipendenti di
% parametro p, uso il seguente ciclo di for:

for i=1:n, x(i)= sum(rand(N,1)<p); end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% VERIFICA EMPIRICA DELLA LEGGE (DEBOLE) DEI GRANDI NUMERI
&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

% La legge debole dei grandi numeri dice che data una successione di
% variabili aleatorie iid X_1, X_2, ... X_n a media finita, la media
% campionaria delle osservazioni converge in probabilità al comune
% valore della media teorica mu = E(X_1). Una verifica empirica di
% questo fatto consiste nel generare n numeri pseudo-aleatorii,
% calcolarne la media campionaria barX= (X_1 + X_2 + ... + X_n)/n quindi
% vedere con un plot come questo valore oscilli intorno a E(X_1) con le
% oscillazioni che diminuiscono al crescere di n.

% Implementazione della verifica:

% 1. Successioni di uniformi su (0,1). In questo caso E(X_1) = 1/2

n=100
x=rand(n,1);
for i=1:n barx(i)=sum(x(1:i))/i;
end

n=1000
x=rand(n,1);
for i=1:n barx(i)=sum(x(1:i))/i;
end

%oppure, anche se concettualmente e` differente

for i=1:100 barx(i) = sum(rand(i,1))/i; end
for i=1:1000 barx(i) = sum(rand(i,1))/i; end
plot(barx,'.')

% Ovviamente, lo stesso si puo' ottenere utilizzando il comando 

n=1000
x=unifrnd(0,1,n,1);
barx(n,1)=0;
barx(1)=x(1)
for i=2:n barx(i)=((i-1)*barx(i-1)+x(i))/i;
end
plot (barx,'.')

% 2. Successioni di bernoulli di parametro(0.3). 
% In questo caso E(X_1) = 0.3


n=1000
x=rand(n,1)<0.3;
barx(n,1)=0;
barx(1)=x(1)
for i=2:n barx(i)=((i-1)*barx(i-1)+x(i))/i;
end
plot (barx,'.')


% Equivalentemente:

n=10000
x=binornd(1,0.3,n,1);
barx(n,1)=0;
barx(1)=x(1)
for i=2:n barx(i)=((i-1)*barx(i-1)+x(i))/i;
end
plot (barx,'.')


% 3. Successioni di binomiali di parametro(200,0.3). In questo caso
% E(X_1) = 200* 0.3 =60

n=1000
x=binornd(200,0.3,n,1);
barx(n,1)=0;
barx(1)=x(1)
for i=2:n barx(i)=((i-1)*barx(i-1)+x(i))/i;
end
plot (barx,'.')

% Provate con altre distribuzioni a vostra conoscenza.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
                  TEOREMA CENTRALE DEL LIMITE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Sia X_1,X_2,........X_n una successione di variabili aleatorie i.i.d. 
% con media  mu e varianza finita sigma^2 e S_n=X_1+X_2+...+X_n la somma.

% Allora, la somma STANDARDIZZATA S_n= (S - n*mu)/(sqrt(n)*sigma)
% converge in distribuzione ad una v.a. normale standardizzata.

% Osservazione.  In generale, la convergenza e' tanto più rapida quanto
% più la distribuzione di partenza e' simmetrica. La regola empirica
% spesso utilizzata e' di applicare l'approssimazione normale se
% n>=30. Ma tale valore deve essere aumentato per leggi fortemente
% asimmetriche, e può essere diminuito nel caso di leggi molto
% simmetriche: ad esempio, per la distribuzione uniforme su un
% intervallo, l'approssimazione è già valida per n>=10.

% Esempio 4: Approssimazione normale della Binomiale(n,p)

% Sia Y~Bin(n,p) con p in [0,1] ed n "grande". Data l'interpretazione di
% Y come somma di n v.a. di bernoulli i.i.d. di parametro p,
% ``asintoticamente'' (che significa per n "sufficientemente" grande)
% Y~N(np, np(1-p)) nel senso che la f.d.r. della variabile
% standardizzata S_n= (Y -np)/sqrt(np(1-p)) calcolata in ogni punto z
% sull'asse reale converge alla f.d.r della legge N(0,1) (ovviamente
% calcolata nello stesso punto La densità binomiale è simmetrica quando
% p=0.5 e tanto più asimmetrica quanto più p è prossimo a 0 od 1. Una
% regola empirica per applicare l'approssimazione normale è la seguente:
% Se np > 5 e n(1-p) > 5 approssima con la legge normale.

% Ricorda  che  con  n ``sufficientemente  grande   e p sufficientemente
% piccolo'', la Binomiale si approssima con la legge Poiss(n*p).


% Osservazione  CORREZIONE DI CONTINUITÀ

% Vale per ogni v.a. discreta. 
% Nel caso della Bin(n,p), si ha:

% Prob(Y <= y) ~ Prob(Z <= (y - np)/sqrt[np*(1-p)]) con Z normale
% standard.  Essendo, Prob(Y <= y) una funzione costante nell'intervallo
% [y,y+1) per y=0,1,2,.....,n (costante "a tratti") mentre Prob(Z <= (y
% - np)/sqrt[np*(1-p)]) è strettamente crescente, si ovvia a questo
% fatto, scegliendo un valore delta in [0,1) (tipicamente delta=0.5) e
% approssimando nel seguente modo: Prob(Y <= y) ~ Prob(Z <= (y + 0.5 -
% np)/sqrt[np*(1-p)] per y=0,1,2,.....n.



% Sovrapponiamo le f.d.r. delle leggi N(0,1) e Bin(n,p), 
% o altre, all'aumentare di n.  

p=0.3;
n=10; 
mu = n*p
sigma = sqrt(n*p*(1-p));
fplot('binocdf(x,10,0.3)',[0,ceil(mu+3*sigma)]);
hold on
% n*p = 3
% sqrt(n*p*(1-p)) = 1.4491
fplot('normcdf(x,3, 1.4491)', [0, mu+3*sigma],'r')
% ATTENZIONE: normcdf vuole la deviazione standard, non la varianza!
hold off

% approssimazione migliore:
p=0.3;
n=30;
mu = n*p
sigma = sqrt(n*p*(1-p));
fplot('binocdf(x,30,0.3)',[0,ceil(mu+3*sigma)]);
hold on 
% n*p = 9
% sqrt(n*p*(1-p)) = 2.51;
fplot('normcdf(x,9, 2.51)', [0, mu+3*sigma],'r')

% Calcolo dell'errore: 
p = 0.05;
nmax = 2000;
i = 1;
for n = 1:10:nmax;
  mu = n*p;
  sigma = sqrt(n*p*(1-p));
  c = ceil (mu+3*sigma);
  b = binocdf (1:c, n, p);
  n1= normcdf (1:c, mu, sigma);       % senza c.d.c.
  n2= normcdf (1.5:c+.5, mu, sigma);  % con c.d.c.
  err1(i) = max (abs (b-n1) );
  err2(i) = max (abs (b-n2) );
  nx(i) = n;
  i = i+1;
end
plot(nx,err1,nx,err2) % visualizza l'andamento dell'errore

% Poissoniana:
n      = 100;
lambda = 50;
x      = [0:n];

y = poisspdf (x, lambda);             
z = normpdf  (x, lambda, sqrt(lambda));    % senza c.d.c.
t = normpdf  (x+.5, lambda, sqrt(lambda)); % con c.d.c.
plot (x, z-y, 'r', x, t-y, 'b');

y = poisscdf (x,lambda);
z = normcdf  (x, lambda, sqrt(lambda));    % senza c.d.c.
t = normcdf  (x+.5, lambda, sqrt(lambda)); % con c.d.c.
plot (x,z-y, 'r',x ,t-y, 'b');


% generiamo una binomiale 
n=5000; % meglio con n=500000, ma matlab e' cosi' lento...
m=15;
p=.4;
for i=1:n, x(i)= sum(rand(m,1)<p); end
z=normrnd(p*m,sqrt(p*m*(1-p)),n,1);
% sovrapponiamole:
nbins = 10;
[n1,x1]=hist(x,nbins);
[n2,x2]=hist(z,nbins);
plot(x1,n1,x2,n2)