###  Guión para el curso de Probabilidad
###  Licenciatura en Matemáticas, UG 
###  Tema 4: cadenas de Markov / Paseo al Azar


#############################################
#############################################
#############################################
#############################################
###  Cadenas de Markov
### Ejemplo del libro de Baclawski modificado
### Paseo al azar con probabilidad de éxito
### 18/38 (ruleta americana), partiendo de 50
### con barreras absorbentes en 0 y 100
### Se grafica la distribución de probabilidad al paso n

p <- 18/38
q <- 1-p
n <- 100
c <- 50
m <- matrix(0,n+1,n+1)
for(i in 2:n){
  m[i,i-1] <- q
  m[i,i+1] <- p
}
m[1,1] <- 1
m[n+1,n+1] <- 1
x <- matrix(0,1,n+1)
x[1,c+1] <- 1
for(i in 1:600){
  x <- x%*% m
  plot((0:n)[x>0], x[x>0], xlim=c(0,n+1),lwd=8,
       ylim=c(0,0.5),xlab='i',ylab='P(X=i)',type='h',
       main='Distribución de probabilidad para un paseo al azar asimétrico')
  abline(v=50,col='red')
  text(25,0.4,paste('n=',i),cex=2)
  Sys.sleep(.1)
}

### n = 1000

p <- 18/38
q <- 1-p
n <- 100
c <- 50
m <- matrix(0,n+1,n+1)
for(i in 2:n){
  m[i,i-1] <- q
  m[i,i+1] <- p
}
m[1,1] <- 1
m[n+1,n+1] <- 1
x <- matrix(0,1,n+1)
x[1,c+1] <- 1
for(i in 1:1000){
  x <- x%*% m
}
plot((0:n)[x>0], x[x>0], xlim=c(0,n+1),lwd=8,
     ylim=c(0,1),xlab='i',ylab='P(X=i)',type='h',
     main='Distribución de probabilidad para un paseo al azar asimétrico')
abline(v=50,col='red')
text(25,0.4,paste('n=',i),cex=2)


### n = 5000
p <- 18/38
q <- 1-p
n <- 100
c <- 50
m <- matrix(0,n+1,n+1)
for(i in 2:n){
  m[i,i-1] <- q
  m[i,i+1] <- p
}
m[1,1] <- 1
m[n+1,n+1] <- 1
x <- matrix(0,1,n+1)
x[1,c+1] <- 1
for(i in 1:5000){
  x <- x%*% m
}
plot((0:n)[x>0], x[x>0], xlim=c(0,n+1),lwd=8,
     ylim=c(0,1),xlab='i',ylab='P(X=i)',type='h',
     main='Distribución de probabilidad para un paseo al azar asimétrico')
abline(v=50,col='red')
text(25,0.4,paste('n=',i),cex=2)







#############################################
###  Trayectorias para el paseo al azar simétrico
set.seed(5678)
nn <- 2000
xi <- 2*rbinom(nn,1,0.5)-1
sn <- cumsum(xi)
plot(1:nn,sn,type='l',ylim=c(-2.5*sqrt(nn),2.5*sqrt(nn)),
     main='Paseo al azar simétrico',xlab='',ylab='')
abline(h=0,col='red')
lines(1:nn, sqrt(2*(1:nn)*log(max(2,log(1:nn)))),
      type='l',col='red')
lines(1:nn, -sqrt(2*(1:nn)*log(max(2,log(1:nn)))),
      type='l',col='red')
for(i in 1:50){
  lines(1:nn,cumsum(2*rbinom(nn,1,0.5)-1),type='l',
        col=rainbow(50)[i])
  Sys.sleep(.1)
}

### reducimos el número de trayectorias
set.seed(5678)
nn <- 2000
xi <- 2*rbinom(nn,1,0.5)-1
sn <- cumsum(xi)
plot(1:nn,sn,type='l',ylim=c(-2.5*sqrt(nn),2.5*sqrt(nn)),
     main='Paseo al azar simétrico',xlab='',ylab='')
abline(h=0,col='red')
ceros <- numeric(5)
ceros[1] <- sum(sn== 0)
lines(1:nn, sqrt(2*(1:nn)*log(max(2,log(1:nn)))),
      type='l',col='red')
lines(1:nn, -sqrt(2*(1:nn)*log(max(2,log(1:nn)))),
      type='l',col='red')
for(i in 1:4){
  sn <- cumsum(2*rbinom(nn,1,0.5)-1)
  lines(1:nn,sn,type='l',col=rainbow(5)[i])
  ceros[i+1] <- sum(sn==0)
  Sys.sleep(.1)
}
legend('topright',legend=ceros,col=c('black',rainbow(5)[1:4]),
       lwd=2)


# Ampliamos las trayectorias a 200,000 juegos
set.seed(5678)
nn <- 200000
xi <- 2*rbinom(nn,1,0.5)-1
sn <- cumsum(xi)
ceros <- numeric(5)
ceros[1] <- sum(sn== 0)
plot(1:nn,sn,type='l',ylim=c(-2.5*sqrt(nn),2.5*sqrt(nn)),
     main='Paseo al azar simétrico',xlab='',ylab='')
abline(h=0,col='red')
lines(1:nn, sqrt(2*(1:nn)*log(max(2,log(1:nn)))),
      type='l',col='red')
lines(1:nn, -sqrt(2*(1:nn)*log(max(2,log(1:nn)))),
      type='l',col='red')
for(i in 1:4){
    sn <- cumsum(2*rbinom(nn,1,0.5)-1)
    lines(1:nn,sn,type='l',col=rainbow(5)[i])
    ceros[i+1] <- sum(sn==0)
    Sys.sleep(.05)
}
legend('topleft',legend=ceros,col=c('black',rainbow(5)[1:4]),
       lwd=2)


### Ruina del jugador
### p = 18/37

###  Graficas de trayectorias
set.seed(65432)
nn <- 10000; pp <- 18/37
sn <- cumsum(c(10,2*rbinom(nn,1,pp)-1))
nstop <- min(c(min((1:nn)[sn==0]),min((1:nn)[sn==20]),nn))
plot(0:(nstop-1), sn[1:nstop],ylim=c(0,20),pch=19,
     xlim=c(0,200),type='l',xlab='',ylab='',main='Paseo al Azar')
points(nstop-1,sn[nstop],pch=19)
abline(h=c(0,20),col='red')
abline(v=0,col='red')

for(i in 1:9){
  sn <- cumsum(c(10,2*rbinom(nn,1,pp)-1))
  nstop <- min(c(min((1:nn)[sn==0]),min((1:nn)[sn==20]),nn))
  lines(0:(nstop-1), sn[1:nstop],type='l',col=rainbow(35+i)[i])
  points(nstop-1,sn[nstop],pch=19,col=rainbow(50)[35+i])
}


### Aprox. probabilidad de ruina a=10, N=20
set.seed(1029)
nncol = 1000; nn=2000; aa=10; NN=20; pp=18/37
funcion0 <- function(x){min(match(0,x),nn,na.rm = TRUE)}
funcionN <- function(x){min(match(NN,x),nn,na.rm = TRUE)}

matt <- matrix(2*rbinom(nn*nncol,1,pp)-1,ncol=nncol)
matt <- apply(matt,2,cumsum)
matt <- rbind(rep(aa,nncol),matt +aa)
vec0 <- apply(matt,2,funcion0)
vecN <- apply(matt,2,funcionN)
sum(vec0<vecN)

### Aprox. probabilidades de ruina para todos los a
set.seed(1029)
nncol = 1000; nn=2000; NN=20; pp=18/37
matt <- matrix(2*rbinom(nn*nncol,1,pp)-1,ncol=nncol)
matt <- apply(matt,2,cumsum)
probruina <- numeric(19)
for(aa in 1:19){
  mattaa <- rbind(rep(aa,nncol),matt +aa)
  vec0 <- apply(mattaa,2,funcion0)
  vecN <- apply(mattaa,2,funcionN)
  probruina[aa] <-  sum(vec0<vecN)/1000
}

plot(1:19,probruina,type='h',lwd=10,xlab='Capital inicial',
     ylab='Probabilidad',ylim=c(0,1),
     main='Probabilidad de Ruina en Ruleta Europea, N=20')

pruina.teo <- function(param,N){
  theta <- (1-param)/param
  pruina <- numeric(N-1)
  for(i in 1:(N-1)){
    pruina[i] <- (theta^N-theta^i)/(theta^N-1)}
  return(pruina)
}    

pruina20 <- pruina.teo(pp,20)
points(1:19,pruina20,pch=19,col='red')
abline(h=(0:10)/10,col='grey')

#############################################
###  Ley del arcoseno. Resultados Teóricos
### Ley discreta

pr0 <- function(k)  # Probabilidad de estar en 0
{  choose(2*k,k)/2^{2*k} }

asd <- function(n,k){ # Ley discreta del arcoseno
  alpha=choose(2*k,k)*choose(2*(n-k),n-k)/2^{2*n}
  return(alpha) }

probs20 <- asd(10,0:9)
cumsum(probs20)

##  Funcion de dist. continua
arcseno <- function(x) {1/(pi*sqrt(x*(1-x)))}
xval <- (1:100)/101
plot(xval,arcseno(xval), ylim=c(0,2.5),type='l',lwd=2,
     xlab='',ylab='',xaxp=c(0,1,10),main='Densidad de la distribución arcoseno')

for(i in 0:10)
  points(i/10,10*asd(10,i),pch=19,col='blue')
for(i in 0:20)
  points(i/20,20*asd(20,i),pch=19,col='red')
legend(0,0.3,c('n=10','n=20'),col=c('blue','red'),
       pch=c(19,19))

### Función de distribución
xval1 <- (0:1000)/1000
plot(xval1,2*asin(sqrt(xval1))/pi,type='l',ylim=c(0,1),lwd=2,
     xlab='',ylab='',main='Función de Distribución del Arcoseno')


#######################################################
###  Ley del arcoseno. Simulación
### Consideramos 20 juegos simétricos

funcion.arc <- function(x){
  pos <- sum(x > 0)
  ceros.arc <- which(x==0)
  pos <-pos + sum(x[ceros.arc-1]>0)
  return(pos)
}
set.seed(7584)
ll <- 20; nncol <- 100000
mat.arc <- matrix(2*rbinom(ll*nncol,1,0.5)-1,ncol=nncol) 
mat.arc <- apply(mat.arc,2,cumsum)
vecpos <- numeric(nn)
vecpos <- apply(mat.arc,2,funcion.arc)
plot(table(vecpos)/nncol,lwd=6, xlab='',ylab='Probabilidad',
     main='Ley del Arcoseno')
points((0:10)*2+.2,asd(10,0:10),col='red',pch=17,cex=1.5)


#######################################################
###  La otra ley del arcoseno
###  Ultima visita al origen

set.seed(192837)
nn = 20; nncol=100000
matt <- matrix(2*rbinom(nn*nncol,1,0.5)-1,ncol=nncol)
matt <- apply(matt,2,cumsum)
matt <- rbind(rep(0,nncol),matt)

funcion.max0 <- function(x){match(0,rev(x))}
vec.arc <- apply(matt,2,funcion.max0)
tab.arc <- (table(vec.arc)-1)/nncol

plot(2*(0:(nn/2)), tab.arc,type='h',ylim=c(0,0.2),
     lwd=4,ylab = '',xlab='',
     main='Ley del arcoseno para la última visita')
points((0:10)*2+.2,asd(10,0:10),col='red',pch=17,cex=1.5)