#####
##### Primera sesion del curso de Probabilidad 2018

#########################################
###   Sucesion de Ensayos de Bernoulli
#########################################


### Gráficas para la distribución Binomial
### Diferentes valores de p, n=20
par(mfrow=c(1,3))
plot(0:20, dbinom(0:20,20,0.1),type='h',lwd=10,col='darkblue', lend=1,
     xlab='Values',ylab='Probabilities', ylim=c(0,.30),
     main='n=20, p=0.1', cex.main=2, cex.lab=1.5)
plot(0:20, dbinom(0:20,20,0.5),type='h',lwd=10,col='darkblue', lend=1,
     xlab='Values',ylab='Probabilities', ylim=c(0,.30),
     main='n=20, p=0.5', cex.main=2, cex.lab=1.5)
plot(0:20, dbinom(0:20,20,0.8),type='h',lwd=10,col='darkblue', lend=1,
     xlab='Values',ylab='Probabilities', ylim=c(0,.30),
     main='n=20, p=0.8', cex.main=2, cex.lab=1.5)
par(mfrow=c(1,1))
### Diferentes valores de n, p=0.5
par(mfrow=c(1,3))
plot(0:10, dbinom(0:10,10,0.5),type='h',lwd=12,col='darkblue', lend=1,
     xlab='Values',ylab='Probabilities', ylim=c(0,.30),main='n=10, p=0.5',
     cex.main = 2, cex.lab=1.5)
plot(0:20, dbinom(0:20,20,0.5),type='h',lwd=10,col='darkblue', lend=1,
     xlab='Values',ylab='Probabilities', ylim=c(0,.30),main='n=20, p=0.5',
     cex.main = 2, cex.lab=1.5)
plot(0:40, dbinom(0:40,40,0.5),type='h',lwd=6,col='darkblue', lend=1,
     xlab='Values',ylab='Probabilities', ylim=c(0,.30),main='n=40, p=0.5',
     cex.main = 2, cex.lab=1.5)
par(mfrow=c(1,1))

##  Poisson, diferentes valores de lambda
par(mfrow=c(1,3))
nn <- 20
plot(0:nn, dpois(0:nn,1),type='h',lwd=10,col='darkblue', lend=1,
     xlab='Values',ylab='Probabilities', ylim=c(0,.4),
     main=expression(paste(lambda,' = 1')),cex.main=2, cex.lab=1.5)
plot(0:nn, dpois(0:nn,5),type='h',lwd=10,col='darkblue', lend=1,
     xlab='Values',ylab='Probabilities', ylim=c(0,.40),
     main=expression(paste(lambda,' = 5')), cex.main = 2, cex.lab=1.5)
plot(0:nn, dpois(0:nn,10),type='h',lwd=10,col='darkblue', lend=1,
     xlab='Values',ylab='Probabilities', ylim=c(0,.40),
     main=expression(paste(lambda,' = 10')), cex.main=2, cex.lab=1.5)
par(mfrow=c(1,1))

##  Geometric (Number of failures before 1st success)
par(mfrow=c(1,3))
nn <- 20
plot(0:nn, dgeom(0:nn,.5),type='h',lwd=10,col='darkblue', lend=1,
     xlab='Values',ylab='Probabilities', ylim=c(0,.5),
     main='p = 0.5',cex.main=2, cex.lab=1.5)
plot(0:nn, dgeom(0:nn,.3),type='h',lwd=10,col='darkblue', lend=1,
     xlab='Values',ylab='Probabilities', ylim=c(0,.50),
     main='p = 0.3', cex.main = 2, cex.lab=1.5)
plot(0:nn, dgeom(0:nn,.1),type='h',lwd=10,col='darkblue', lend=1,
     xlab='Values',ylab='Probabilities', ylim=c(0,.50),
     main='0 = 0.1', cex.main=2, cex.lab=1.5)
par(mfrow=c(1,1))

##  Geometric (Number of trials until 1st success)
par(mfrow=c(1,3))
nn <- 20
plot(0:nn+1, dgeom(0:nn,.5),type='h',lwd=10,col='darkblue', lend=1,
     xlab='Values',ylab='Probabilities', ylim=c(0,.5),xlim=c(0,20),
     main='p = 0.5',cex.main=2, cex.lab=1.5)
plot(0:nn+1, dgeom(0:nn,.3),type='h',lwd=10,col='darkblue', lend=1,
     xlab='Values',ylab='Probabilities', ylim=c(0,.50),xlim=c(0,20),
     main='p = 0.3', cex.main = 2, cex.lab=1.5)
plot(0:nn+1, dgeom(0:nn,.1),type='h',lwd=10,col='darkblue', lend=1,
     xlab='Values',ylab='Probabilities', ylim=c(0,.50),xlim=c(0,20),
     main='0 = 0.1', cex.main=2, cex.lab=1.5)
par(mfrow=c(1,1))

# Negative Binomial (Number of failures before kth success)
par(mfrow=c(1,3))
plot(0:40, dnbinom(0:40,2,0.5),type='h',lwd=6,col='darkblue', lend=1,
     xlab='Values',ylab='Probabilities', ylim=c(0,.30),
     main='k=2, p=0.5', cex.main=2, cex.lab=1.5)
plot(0:40, dnbinom(0:40,5,0.5),type='h',lwd=6,col='darkblue', lend=1,
     xlab='Values',ylab='Probabilities', ylim=c(0,.30),
     main='k=5, p=0.5', cex.main=2, cex.lab=1.5)
plot(0:40, dnbinom(0:40,10,0.5),type='h',lwd=6,col='darkblue', lend=1,
     xlab='Values',ylab='Probabilities', ylim=c(0,.30),
     main='k=10, p=0.5', cex.main=2, cex.lab=1.5)
par(mfrow=c(1,1))

# Negative Binomial (Number of trials until kth success)
par(mfrow=c(1,3))
plot(0:40+2, dnbinom(0:40,2,0.5),type='h',lwd=5,col='darkblue', lend=1,
     xlab='Values',ylab='Probabilities', ylim=c(0,.30),xlim=c(0,50),
     main='k=2, p=0.5', cex.main=2, cex.lab=1.5)
plot(0:40+5, dnbinom(0:40,5,0.5),type='h',lwd=5,col='darkblue', lend=1,
     xlab='Values',ylab='Probabilities', ylim=c(0,.30),xlim=c(0,50),
     main='k=5, p=0.5', cex.main=2, cex.lab=1.5)
plot(0:40+10, dnbinom(0:40,10,0.5),type='h',lwd=5,col='darkblue', lend=1,
     xlab='Values',ylab='Probabilities', ylim=c(0,.30),xlim=c(0,50),
     main='k=10, p=0.5', cex.main=2, cex.lab=1.5)
par(mfrow=c(1,1))



set.seed(192837)
index <- 1:100

par(mfrow=c(4,1))
plot(index,rbinom(100,1,0.1), type='h',lwd=3,
     main='Bernoulli Trials p=0.1',ylab='',
     bty='n',yaxt='n')
plot(index,rbinom(100,1,0.3), type='h',lwd=3, main='Bernoulli Trials p=0.3',
     ylab='',bty='n',yaxt='n')
plot(index,rbinom(100,1,0.5), type='h',lwd=3, main='Bernoulli Trials p=0.5',
     ylab='',bty='n',yaxt='n')
plot(index,rbinom(100,1,0.7), type='h',lwd=3, main='Bernoulli Trials p=0.7',
     ylab='',bty='n',yaxt='n')
par(mfrow=c(1,1))

### Distribución binomial 
### La función rbinom(k,n,p) simula k valores de una binomial
### con parámetros n y p

rbinom(10, 20, 0.2)
nn <- 1000
vals <- tabulate(rbinom(nn,20,0.2)+1,21)
plot(0:20, vals/nn,pch=16, cex.main=1.7, type='h', lwd=8,
     main='Binomial probability function with n=20, p = 0.2',
     xlab='Values',ylab='Frequency',cex.lab=1.5)
abline(h=0, col='red')
points(0:20, dbinom(0:20,20,0.2), pch=16,col='red')


###  Distribución geométrica
muestra <-  rbinom(100,1,0.1)
index[muestra==1][1] # Ubicacion del primer exito
index[muestra==1][3] 

# Construimos una funcion para hacer esto
geo1 <-  function(n=1000,p=0.1)
  (1:n)[rbinom(n,1,p)==1][1]

geo1()
replicate(10,geo1())

values <- tabulate(replicate(10000,geo1()),50)
plot(1:50,values/10000,pch=16, cex.main=1.7,type='h',lwd=4,
     main='Geometric probability function with p = 0.1',
     xlab='Values',ylab='Frequency',cex.lab=1.5)
abline(h=0, col='red')
abline(v=0, col='red')

points(1:50,dgeom(0:49,0.1),col='red',pch=18)
legend('topright',c('Empirical frequency','Theoretical probability'),
       pch=c(16,18),col=c('black','red'))

###  En R la funcion rgeom genera muestras de la distribucion geometrica
plot(0:50,tabulate(rgeom(10000,0.1)+1,51)/10000,pch=16,ylab='Frequency',
     main='Geometric probability function with p = 0.1',
     xlab='Value',cex.lab=1.5,cex.main=1.7,ylim=c(0,0.1))
abline(h=0, col='red')
abline(v=0, col='red')
points(0:50,dgeom(0:50,0.1),col='red',pch=18)
legend('topright',c('Empirical frequency','Theoretical probability'),
       pch=c(16,18),col=c('black','red'))

### pero no es la misma distribucion geometrica que simulamos antes
### Esta cuenta el numero de fracasos hasta el primer exito
### SIN INCLUIR el éxito.
### Para que coincida con distribucion anterior tenemos que correr los 
### valores a la derecha una unidad
plot(1:50,tabulate(rgeom(10000,0.1),50),pch=16,xlab='k',ylab='Frecuencia',
     main='10,000 ensayos de la distribución geométrica')
abline(h=0, col='red')
abline(v=0, col='red')

xval <- 0:40; clrs <- rainbow(50)
plot(xval+1, dgeom(xval,0.8), type='l',main='Densidad de la Distribución Geométrica',lwd=2,
     ylab='Densidad', xlab='x',col=clrs[4])
lines(xval+1, dgeom(xval,0.5),col=clrs[14],lwd=2)
lines(xval+1, dgeom(xval,0.3),col=clrs[34],lwd=2)
lines(xval+1, dgeom(xval,0.1),col=clrs[44],lwd=2)
legend('topright', c('0.7','0.5','0.3','0.1'),lwd=rep(2,4),
       col=c(clrs[4],clrs[14],clrs[34],clrs[44]))
abline(h=0, col='red')
abline(v=0, col='red')

## Los colores los escogimos con rainbow:
pie(rep(1, 100), col = rainbow(100))

op <- par(no.readonly = TRUE) # the whole list of settable par's.
par(mfrow=c(2,2))
par(mar=c(4,4,4,2)+0.1)
plot(xval+1, dgeom(xval,0.7), type='h',main='Geometric probability function, p = 0.7',lwd=3,
     ylab='Probability', xlab='x', col = clrs[36])
plot(xval+1, dgeom(xval,0.5),col=clrs[34],lwd=3,type='h',
     main='Geometric probability function, p = 0.5',ylab='Probability', xlab='x')
par(mar=c(5,4,3,2)+0.1)
plot(xval+1, dgeom(xval,0.3),col=clrs[32],lwd=3,type='h',
     main='Geometric probability function, p = 0.3',ylab='Probability', xlab='x')
plot(xval+1, dgeom(xval,0.1),col=clrs[30],lwd=3,type='h',
     main='Geometric probability function, p = 0.1',ylab='Probability', xlab='x')
par(mfrow=c(1,1))
par(op)


###  Binomial negativa
###  Podemos construirla como suma de k geometricas
set.seed(1234)
sum(replicate(3,geo1()))

bineg <-  function(k=2,nn=1000,pp=0.1)
  (1:nn)[rbinom(nn,1,pp)==1][k]
bineg()
bineg(10)

valuesnb <- tabulate(replicate(10000,bineg(3)),100)
plot(1:100,valuesnb,pch=16,main='Negative Binomial with p = 0.1 and k = 3',
     xlab='Values',ylab='Frequency', cex.main=1.5, cex.lab=1.5)
abline(h=0, col='red')
abline(v=0, col='red')

##  Densidades binomial negativa


op <- par(no.readonly = TRUE) # the whole list of settable par's.
par(mfrow=c(2,2))
par(mar=c(4,4,4,2)+0.1)
plot(xval+3, dnbinom(xval,3,0.7), type='h',main='Distribución Binomial Negativa, k=3, p=0.7',lwd=2,
     ylab='Densidad', xlab='x',xlim=c(0,50))
plot(xval+3, dnbinom(xval,3,0.5),col=clrs[4],lwd=2,type='h'
     ,main='Distribución Binomial Negativa, k=3, p=0.5',ylab='Densidad', xlab='x',xlim=c(0,50))
par(mar=c(5,4,3,2)+0.1)
plot(xval+3, dnbinom(xval,3,0.3),col=clrs[28],lwd=2,type='h',
     main='Distribución Binomial Negativa, k=3, p=0.3',ylab='Densidad', xlab='x',xlim=c(0,50))
plot(xval+3, dnbinom(xval,3,0.1),col=clrs[34],lwd=2,type='h',
     main='Distribución  Binomial Negativa, k=3, p=0.1',ylab='Densidad', xlab='x',xlim=c(0,50))

par(mar=c(4,4,4,2)+0.1)
plot(xval+3, dnbinom(xval,3,0.3), type='h',main='Distribución Binomial Negativa, k=3, p=0.3',lwd=2,
     ylab='Densidad', xlab='x',xlim=c(0,50))
plot(xval+5, dnbinom(xval,5,0.3),col=clrs[4],lwd=2,type='h'
     ,main='Distribución Binomial Negativa, k=5, p=0.3',ylab='Densidad', xlab='x',xlim=c(0,50))
par(mar=c(5,4,3,2)+0.1)
plot(xval+7, dnbinom(xval,7,0.3),col=clrs[28],lwd=2,type='h',
     main='Distribución Binomial Negativa, k=7, p=0.3',ylab='Densidad', xlab='x',xlim=c(0,50))
plot(xval+9, dnbinom(xval,9,0.3),col=clrs[34],lwd=2,type='h',
     main='Distribución  Binomial Negativa, k=9, p=0.3',ylab='Densidad', xlab='x',xlim=c(0,50))

par(mfrow=c(1,1))
par(op)

################################
###  Distribución binomial
################################

# Grafica para n=10 y p variando de 0 a 1 en 0.01
binom <- function(p){
  plot(0:10,dbinom(0:10,10,p),type='h',lwd=5,ylim=c(0,0.5),
       xlab='Valores',ylab='Probabilidad', main=p)
  Sys.sleep(0.1)
}
ignore <- sapply((0:100)/100,binom)

# Gráfica para p=1/2 y n variando de 1 a 100

binom2<- function(n){
  plot(0:n,dbinom(0:n,n,0.5),type='h',lwd=5,ylim=c(0,0.5),
       xlab='Valores',ylab='Probabilidad',main=n)
  Sys.sleep(0.1)
}
ignorar <- sapply((0:100),binom2)

## Cambiamos la escala
binom3<- function(n){
  plot(0:n,dbinom(0:n,n,0.5),type='h',lwd=5,ylim=c(0,0.5),
       xlab='Valores',ylab='Probabilidad',main=n,
       xlim=c((n/2)-2*sqrt(n),(n/2)+2*sqrt(n)))
  Sys.sleep(0.08)
}
ignorar <- sapply((0:200),binom3)

## Cambiamos p
binom4<- function(n){
  plot(0:n,dbinom(0:n,n,0.25),type='h',lwd=5,ylim=c(0,0.5),
       xlab='Valores',ylab='Probabilidad',main=n,
       xlim=c((n/4)-2*sqrt(n),(n/4)+2*sqrt(n)))
  Sys.sleep(0.08)
}
ignorar <- sapply((0:200),binom4)