# 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)-mu^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)