# funzioni di probabilita'
 
# poissoniana
dpois(0,1)
exp(-1)
 
# binomiale  --------------------------------------------------------
dbinom(5, 10, 0.5)
dbinom(500000, 1000000, 0.5)
 
n<-10; p<-0.2
dbinom(0:n, n, p)
pbinom(0:n, n, p)
 
# grafico
n=10; p=0.5; plot( 0:n, dbinom(0:n, n, p), ty='h' )
n=10; p=0.2; plot( 0:n, dbinom(0:n, n, p), ty='h' )
 
# meglio ancora
n=10; p=0.5; barplot(dbinom(0:n, n, p), names=c(0:n), col='cyan', xlab='x', ylab='f(x)', main=sprintf('Binomiale n=%.0f, p=%.2f',n,p))
 
# plot a gradini della cumulativa:
n=10; p=0.2; plot( 0:n, pbinom(0:n, n, p), ty='s' )
 
# calcolo media e deviazione standard:
n=10; p=0.2; x=c(0:n); fx=dbinom(x, n, p); mu=sum(x*fx); sigma=sqrt(sum(x^2*fx)-m^2)
mu
sigma
# confrontiamo con
n*p
sqrt(n*p*(1-p))
 
# Geometrica -----------------------------------------------------------
p=1/8; pgeom(1,p)
# confrontiamo con
p
# Att: la geometrica di R e' definita non come tentativo
# al quale si ha il primo successo, ma nr di tentativi _prima_di_ 
# avere un successo, quindi comincia da 0.
 
 
# Gaussiana  ----------------------------------------------------------
dnorm(0)             # la gaussiana ha default mu=0, sigma=1
pnorm(0)
pnorm(1)-pnorm(-1)
 
pnorm(7,5,2)-pnorm(3,5,2)
 
1-pnorm(3)
pnorm(Inf)-pnorm(3)
 
 
qnorm(0.25)          # quantili
qnorm(1)
qnorm(0)
qnorm(0.75)-qnorm(0.25)
qnorm( seq(0, 1, by=0.1) ) 
 
# grafico di pdf e cumulativa
mu=10; sigma=3
fx <- function(x) dnorm(x,mu,sigma)
Fx <- function(x) pnorm(x,mu,sigma)
layout( matrix(c(1:2),2, 1) )
plot(fx, mu-4*sigma, mu+4*sigma, ty='l', col='cyan',
     main=sprintf('Gaussiana con mu=%.1f, sigma=%.1f [pdf]', mu, sigma))
plot(Fx, mu-4*sigma, mu+4*sigma, ty='l', col='blue',
     main=sprintf('Gaussiana con mu=%.1f, sigma=%.1f  [cumulativa]', mu, sigma))
layout(1)
 
# campione estratti secondo una gaussiana 
x<-rnorm(1000)
length(x[x>=0])
length(x[x>=-1 & x<=1])
length(x[abs(x)<=1])
 
mean(x)
sd(x)
 
hist(x, nc=50, prob=TRUE) 
plot(dnorm, -4, 4, ty='l', add=TRUE, col='red')
fx <- function(z) dnorm(z, mean(x), sd(x))     
plot(fx, -4, 4, ty='l', add=TRUE, col='brown')
 
# facciamoci 'in casa' media e deviazione standard:
my.mean <- function(x) sum(x)/length(x)
my.sd   <- function(x) sqrt( my.mean(x^2) - my.mean(x)^2 )
my.mean(x) #  = mean(x), OK
my.sd(x)   # != sd(x)  , ?!
# indaghiamo:
my.sd(x)*sqrt(length(x)/(length(x)-1)) 
# sd() e' definita in questo modo, per motivi cari agli
# statistici: INTERFERENZA fra statistica descrittiva e statistica
# inferenziale: nonsense!
 
 
# somma di variabili gaussiane ------------------------------------
x1<-rnorm(10000)
x2<-rnorm(10000)
x3<-rnorm(10000)
x4<-rnorm(10000)
x5<-rnorm(10000)
 
xt<-x1^2+x2^2+x3^2+x4^2+x5^2
mean(xt)
sd(xt)
 
c2<-rchisq(10000,5)
mean(c2)
sd(c2)
 
#-- ancora sulla geometrica ---------------------------------------
dgeom(1,1/18)
dgeom(18,1/18)
pgeom(18,1/18)
 
pr180 <- 1-pgeom(180,1/18)
# Immaginiamo di pensare alle PROSSIME 180 estrazioni
# quanto vale la probabilita' che un numero non esca?
# (90 numeri * 10 ruote):
dbinom(0, 900, pr180)
dbinom(1, 900, pr180)   # 2.8% : niente di sorprendente! 
dbinom(2, 900, pr180)