###  Guión para el Tema 2: Ensayos de Bernoulli
###  Curso de Probabilidad
###  Universidad de Guanajuato

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

set.seed(192837)
index <- 1:100
muestra <- rbinom(100,1,0.1)

plot(index,muestra, type='h',lwd=2)

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


###  Distribución geométrica
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())

valores <- tabulate(replicate(10000,geo1()),50)
plot(1:50,valores,pch=16, main='Distribución Geométrica con p=0.1')
abline(h=0, col='red')
abline(v=0, col='red')


###  En R la funcion rgeom genera muestras de la distribucion geometrica
plot(0:49,tabulate(rgeom(10000,0.1),50),pch=16,xlab='k',ylab='Frecuencia',
     main='Distribución Geométrica con p=0.1')
abline(h=0, col='red')
abline(v=0, col='red')
### pero no es la misma distribucion geometrica que simulamos antes
### Esta cuenta el numero de fracasos hasta el primer exito
### SIN INCLUIR el exito.
### 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, 50), col = rainbow(50))

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='Distribución Geométrica, p=0.7',lwd=2,
     ylab='Densidad', xlab='x')
plot(xval+1, dgeom(xval,0.5),col=clrs[4],lwd=2,type='h'
     ,main='Distribución Geométrica, p=0.5',ylab='Densidad', xlab='x')
par(mar=c(5,4,3,2)+0.1)
plot(xval+1, dgeom(xval,0.3),col=clrs[28],lwd=2,type='h',
     main='Distribución Geométrica, p=0.3',ylab='Densidad', xlab='x')
plot(xval+1, dgeom(xval,0.1),col=clrs[34],lwd=2,type='h',
     main='Distribución Geométrica, p=0.1',ylab='Densidad', 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)

valoresbn <- tabulate(replicate(10000,bineg(3)),100)
plot(1:100,valoresbn,pch=16,main='Binomial negativa com p=0.1 y k=3',
     xlab='k',ylab='Frecuencia')
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')
  Sys.sleep(0.1)
}
ignorar <- 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')
  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',
       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',
       xlim=c((n/4)-2*sqrt(n),(n/4)+2*sqrt(n)))
  Sys.sleep(0.08)
}
ignorar <- sapply((0:200),binom4)