% Terzo  laboratorio 29 ottobre 2002
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Argomento: variabili aleatorie discrete

% Impariamo ad usare matlab per calcolare funzione di densità (pdf:
% probability density function), funzione di ripartizione (cdf:
% cumulative distribution function) e quantili di variabili aleatorie
% discrete "notevoli": Binomiale e Poisson.

% Esercizio 1: Sia X una variabile aleatoria discreta con distribuzione
% di Poisson di parametro 3. 
% Calcolare 1) P(X=3), 2) P(X>=3), 3) P(1<=X<5), 4) P(X>=15)

% Sintassi: 
% densità di probabilità   della Poisson(lambda): P(X=x)  = poisspdf(x,lambda)
% funzione di ripartizione della Poisson(lambda): P(X<=x) = poisscdf(x,lambda) 

% Risp 1:    
p1 = poisspdf(3,3);                 % p1 = 0.2240
p2 = 1- poisscdf(2,3);              % p2 = 0.5768
p3 = poisscdf(4,3) - poisscdf(0,3); % p3 = 0.7655
p4 = 1- poisscdf(14,3);             % p4 = 6.7039e-07


% Esercizio 2 Ripetere l'esercizio 1 nel caso in cui X abbia legge
% binomiale di parametri (n=18,p=1/3)

% Sintassi: 
% densità di probabilità   della B(N,p): P(X=x)  = binopdf(x,N,p)
% funzione di ripartizione della B(N,p): P(X<=x) = binocdf(x,N,p) 

% Risp 2 
p1 = binopdf(3,18,1/3);
p2 = 1- binocdf(2,18,1/3);
p3 = binocdf(4,18,1/3) - binocdf(0,18,1/3);
p4 = 1- binocdf(14,18,1/3);


% Esercizio 3.  All'inizio di ogni mese un negoziante fa gli ordini di
% un certo prodotto per le vendite di tutto il mese. Egli può
% plausibimente assumere che il numero delle unità di prodotto vendute
% durante il mese segua una legge di Poisson di parametro lambda=4.

% a. Quanto vale la probabilità che venda al massimo tre unità del
% prodotto in un mese?

% Risp: 

p3 = poisscdf (3, 4);

% b.  Qual è il numero minimo di unità da ordinare perche' la
% probabilità di non rimanere senza prodotto sia almeno 0.9?

% Risp: Dobbiamo determinare q tale che: P(X<=q) >=0.9

% In altri termini, dobbiamo calcolare il quantile di ordine 0.9 della
% legge di Poisson di parametro 4:

% Sintassi:

qalpha = poissinv(alpha,LAMBDA) % alpha è l'ordine del quantile che
                                % stiamo cercando

x = poissinv(0.9,4); % x = 7


% Esercizio 4. Disegnare il grafico della densità di probabilità della
% legge binomiale  di parametri n=40 e p = 7/10, nonchè della f. di ripartizione 

% non funziona perche' pdf ha valori nulli per x non intero
% fplot('binopdf(x,40,7/10)',[0,40])
% fplot('binocdf(x,40,7/10)',[0,40],'o')

% oppure (meglio)
x = 0:40
y = binopdf(x,40,7/10)
bar(x,y)
plot(x,y,'o-')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Calcolo valori di sintesi: Per il calcolo di media e varianza,
% concentriamo l'attenzione sulla binomiale B(n,p).

% il comando è 
[mu,sigma] = binostat(n,p)

% Possiamo calcolare simultaneamente le coppie (mu,sigma) al variare di
% n: Per esempio: con il comando n = linspace(10,100,N) 
% generiamo un insieme di N punti da 10 a 100 egualmente spaziati.

% n = linspace(10,100,10)
% oppure:
n = [10:10:100]

% n =
% 
% 10    20    30    40    50    60    70    80    90   100

% A questo punto, con 
[mu,sigma] = binostat(n,0.5)

% fornisce
% mu =  5    10    15    20    25    30    35    40    45    50
% sigma  =  2.5000    5.0000    7.5000   10.0000   12.5000   15.0000   17.5000
      % 20.0000   22.5000   25.0000

% Verifica empirica del seguente risultato:

% Date due variabili aleatorie X, Y con funzione di ripartizione F, G
% rispettivamente, allora E(X) < E(Y) se e solo se F(x)>G(x) per ogni x
% scelto fra i reali. Esempio:

fplot('binocdf(x,40,7/10)',[0,40])
hold on 
fplot('binocdf(x,40,8/10)',[0,40])


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% VERIFICA EMPIRICA DI COME POISSON APPROSSIMI BINOMIALE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% criterio empirico: n/p > 500
% Sia lambda = 5

% 1) 5 = 20 * 1/4    n/p = 80
% 2) 5 = 30 * 1/6    n/p = 180
% 3) 5 = 50 * 1/8    n/p = 320
% 4) 5 = 100* 1/20   n/p = 2000

x = 0:20
z = poisspdf(x,5)
plot(x,z,'k')

hold on
y1 = binopdf(x,20,1/4)
plot (x, y1,'b')

y2 = binopdf(x,30,1/6)
plot(x,y2,'y')

y3 = binopdf(x,40,1/8)
plot (x,y3,'g')

y4 = binopdf(x,100,1/20)
plot(x,y4,'r')

hold off

plot(x,z-y1,x,z-y2,x,z-y3,x,z-y4)

max(abs(z-y1))
max(abs(z-y2))
max(abs(z-y3))
max(abs(z-y4))

%-------------------------------------------------------------------------------
% Geometrica

% supponiamo che la probabilita' che un dato componente
% elettronico cessi di funzionare in un dato giorno e' p=0.002.
% Qual'e' la probablita' che esso duri esattamente 50 giorni?

p = 0.002;
d = 50;
prob = geopdf (d, p) % 0.0018

% Che duri 500 giorni o meno?
% sum_{i=0}^{500} geopdf (i, p)
ngiorni = 500;
prob = 0;
for i = 0:ngiorni
   prob = prob+geopdf(i,p);
end
prob
% oppure:
prob = geocdf (ngiorni, p)

% dopo quanti giorni la probabilita' che si sia guastato in quel 
% lasso di tempo vale 1/2?
geoinv (0.5, p) % 346

% calcolo di media e varianza:
[m,v] = geostat(p) % m=499 v=249500 dev_std=499.5

% plot della f.d.r.:
xmax = 1000
x = [0:xmax];
y = geocdf (x,p);
plot (x, y)
% oppure:
fplot ('geocdf(x,.002)', [0 xmax]);
% o anche:
plot ([0:xmax], geocdf([0:xmax],p));
%-------------------------------------------------------------------------------


% per la normale di parametri mu, sigma:
P = normpdf(x,mu,sigma)
P = normcdf(x,mu,sigma) 
alpha=0.4
q = norminv(alpha,mu,sigma)      

% per l'esponenziale di parametro lambda:
P = exppdf(x,1/lambda) 
P = expcdf(x,1/lambda)
alpha=0.4
q = expinv(alpha,1/lambda) % (quantile di ordine alpha)  

% per la uniforme continua sull'intervallo (a,b)
a=1;
b=10;
P= unifpdf(x,a,b)
P= unifcdf(x,a,b)
alpha=0.4
q= unifinv(alpha,a,b)