“Paradoja” del cumpleaños

Seguramente haz escuchado de la “paradoja” del cumpleaños donde se desea determinar la probabilidad de que dos personas en un salón cumplan el mismo día. Para fines de este ejercicio considera que se \(n\) personas, los años bisiestos no son contados ni admiten las personas gemelas; además de que los posibles \(365\) cumpleaños tienen la misma probabilidad de ocurrir.

En resumen, se tienen las siguentes expresiones para determinar la probabilidad, bajo las condiciones anteriores, de que dos personas cumplan el mismo día y de que otra persona cumpla el mismo día que tu.

\[\begin{array}{ccc} \mathbb{P}= \left\{\begin{array}{ll}1-\frac{365!}{365^n(365-n)!} & 1\leq n \leq365\\ 1 & n>365 \end{array}\right. &; &\mathbb{P} = 1-\left(\frac{364}{365}\right)^n \end{array}\]
  1. Crea una función, que de acuerdo a una \(n\) válida, determine ambas probabilidades.
  2. Crea una gráfica donde se tengan la distribución de cada una de las probabilidades y determinar si existe algún momento en el que hay la misma probabilidad, para una \(n\), de que dos personas cumplan el mismo día y de que otra persona cumpla el mismo día que tu.

Podemos simplificar la primera probabilidad como \(\mathbb{P}=1-\frac{365!}{365^n \left(365-n\right)} =1-\left(\frac{365}{365} \ldots \frac{365-n+1}{365}\right)\). Esto con el fin de desasernos de \(365!\) pues es un número tan grande que R lo toma como Inf y no podemos hacer operaciones con ello.

Comencemos por calculara las probabilidades Para la primera que es Para un grupo de n personas la probabilidad de que dos personas cumplan años el mismo día

n<- numeric
prob1 <- function(n){
 x <- 1
 for (i in 1:n){
   if(i>=1 & i<=365){
  
     x = (x*(365-(i-1))/365)
    p1=1-x
     
  } else if(i>365){
    p1<-1
  } else if(i<1){
    p1<-0
  }
 }

 return(p1)
}

prob1(80)
## [1] 0.9999143

Para la segunda probabilidad tenemos Para un grupo de n personas la probabilidad que alguien cumpla años el mismo día que tu

n <- numeric
prob2 <- function(n){
  p2<- numeric
  p2<- 1 -(364/365)^n
  return(p2)
}
prob2(80)
## [1] 0.1970629

Para grficar, tenemos que hacer

library(ggplot2)
vector1<-c()
vector2<-c()
for(i in 1:365){
  vector1[i]<-prob1(i)
  vector2[i]<-prob2(i)
}

data <- data.frame(vector1, vector2, prob2(1:365))

ggplot(data, aes(1:365))+ geom_point(aes(y=vector1), col="red") + geom_point(aes(y=vector2), col="green")

Vamos aumentar la n, a ver que pasa

vector1<-c()
vector2<-c()
for(i in 1:2000){
  vector1[i]<-prob1(i)
  vector2[i]<-prob2(i)
}

data <- data.frame(vector1, vector2, prob2(1:2000))

ggplot(data, aes(1:2000))+ geom_line(aes(y=vector1), col="red") + geom_line(aes(y=vector2), col="green")

Vemos que las probabilidades son casi iguales conforme la \(n\) se acerca a cero y tambien son muy parecidos si hacemos que \(n\) se vaya a infinito.

Relación Fibonnaci-Eigen (Vectores/Valores)

Existen aplicaciones muy interesantes donde se utilizan los conocidos eigen/valores de una matriz.Una de ellas es la relación que tienen estos con los conocidos números de Fibonnacci. Recuerda que los números de Fibonacci quedan representados por la ecuación recursiva \(F_{n} = F_{n-1} + F_{n-2}\) y de una manera muy sencilla se pude ver que

\[\begin{pmatrix} F_n \\F_{n-1} \end{pmatrix} = \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix} \begin{pmatrix} F_{n-1}\\ F_{n_2} \end{pmatrix}\]

  1. Crea una función para obtener el \(n\)-ésimo número de Fibonacci.
  2. Determina mediante el uso de R el eigen valor positivo correspondiente a dicha matriz (es decir, el famoso número áurio o número de oro).
  3. Crea una gráfica para un \(n\) que desees, donde cada punto corrresponda a \((F_{n−1}, ~ F_{n−2})\) o \((F_{n}, ~ F_{n−1})\). Dichos puntos deben ser de color negro
  4. En la misma gráfica coloca la recta que pasa por el eigen vector correspondinete al eigen valor del punto 1.
  5. Elige algún punto de los graficados en el punto dos y multiplicalo por el eigen valor del punto 1 y graficalo en color rojo ¿Qué sucedio?
  6. ¿Qúe concñuyes de todo esto? Veamos que la sucesion de Fibonacci comienza con el \(0\) y el \(1\), esto es que \(F_0=0\) y \(F_1=1\) entonces apartir de \(n=2\) vamos a usar la formula \(F_n=F_{n-1}+F_{n-2}\) y para encontrarlo vamos a usar el siguiente codigo
suc_fibo <- function(n){
  if(n==0){
    return(0)
  } else if(n==1){
    return(n)
  }
  else if(n>=2){
    return(suc_fibo(n-1)+suc_fibo(n-2))
  } else if(n<0){
    cat("No es posible calcular para", n, "intenta con un número mayor o igual a cero")
  }
}

Vamos a probarla

suc_fibo(0)
## [1] 0
suc_fibo(1)
## [1] 1
suc_fibo(5)
## [1] 5
suc_fibo(10)
## [1] 55
suc_fibo(-5)
## No es posible calcular para -5 intenta con un número mayor o igual a cero

Ahora vamos a determinar el número áureo o número de oro, que esta relacionado con la suecesión de Fibonacci, ya si hagarramos una \(n\) cada vez más grande \(\frac{F_{n+1}}{F_{n}}\) se va acercando al número áureo. Para calcular el vaor dada una matriz, vamos a usar

X<-matrix(
  c(
    1, 1,    
    1, 0
  ),
  nrow = 2, byrow = TRUE   
)

#Usamos la función eigen para obtener el número aureo
eigen(X)$values
## [1]  1.618034 -0.618034

Vemos que el primer valor es el que se parce al número buscado entonces

oro <- eigen(X)$values[1]
oro
## [1] 1.618034

Ahora graficaremos los primeros 20 números de la secesión de Fibonacci, es decir \(n=10\)

library(ggplot2)
vec1<-c()
vec2<-c()
for(i in 1:10){
  vec1[i+1]<-suc_fibo(i)
  vec2[i+1]<-suc_fibo(i-1)
}

data1 <- data.frame(vec1,vec2)



ggplot(data1, aes(vec1))+geom_point(aes(y=vec2), col="black")
## Warning: Removed 1 rows containing missing values (geom_point).

Le agregamos una recta que pase por el eigen valor encontrado

library(ggplot2)
vec1<-c()
vec2<-c()
for(i in 1:10){
  vec1[i+1]<-suc_fibo(i)
  vec2[i+1]<-suc_fibo(i-1)
}

data1 <- data.frame(vec1,vec2)

eigen(X)$vector
##            [,1]       [,2]
## [1,] -0.8506508  0.5257311
## [2,] -0.5257311 -0.8506508
ggplot(data1, aes(vec1))+geom_point(aes(y=vec2), col="black")+geom_abline(slope =eigen(X)$vector[2,1]/eigen(X)$vector[1,1])
## Warning: Removed 1 rows containing missing values (geom_point).

Ahora agreguemos el punto

library(ggplot2)
vec1<-c()
vec2<-c()
for(i in 1:10){
  vec1[i+1]<-suc_fibo(i)
  vec2[i+1]<-suc_fibo(i-1)
}

data1 <- data.frame(vec1,vec2)
data1
##    vec1 vec2
## 1    NA   NA
## 2     1    0
## 3     1    1
## 4     2    1
## 5     3    2
## 6     5    3
## 7     8    5
## 8    13    8
## 9    21   13
## 10   34   21
## 11   55   34
eigen(X)$vecto
##            [,1]       [,2]
## [1,] -0.8506508  0.5257311
## [2,] -0.5257311 -0.8506508
ggplot(data1, aes(vec1))+geom_point(aes(y=vec2), col="black")+geom_abline(slope =eigen(X)$vector[2,1]/eigen(X)$vector[1,1]) +geom_point(aes(x=data1[5,1]*oro, y=data1[5,2]*oro),col="red")
## Warning: Removed 1 rows containing missing values (geom_point).

Esta es una forma de ver que si damos un punto en la sucesion de Fibonacci ya la muliplicamos por el valor del numero de oro casi obtenemos otro valor de la serie de Fibonacci, por lo que si aumentamos la n, la misma sucesion de fibonacci nos podrá decir el número áureo.

Procesamiento de textos

La idea será sencilla aunque la implementación no lo sé, por lo que podrías obtener más resultados de los que coloco aquí. De acuerdo a un archivo .txt que se te será proporcionado determinar lo siguiente:

  1. Cantidad de letras.
  2. Cantidad de vocales. 3.Cantidad de espacios. 4.Porcentaje que representa cada letra en el texto. No importa si consideras a los caracteres especiales como letras o no, tampoco si haces distinción entre mayúsculas y minúsculas, ni tampoco las veces que tengas que cargar el archivo.

Vamos a analizar un texto de Romeo y Julieta, para ello descargamos la paqueteria de \(stringr\)

library(stringr)
texto <- readr :: read_file("Romeo_and_Juliet.txt")
let <- sum(str_count(texto,letters))
voc <- str_count(texto, "a|e|i|o|u|A|E|I|O|U")
espa <- str_count(texto, " ")

porcen <-c()
for(i in 1:length(letters)){
  porcen[i] <- str_count(texto, letters[i])/let*100
}
cat("En el texto hay", let, "letras")
## En el texto hay 95140 letras
cat("En el texto hay", voc, "vocales")
## En el texto hay 40531 vocales
cat("En el texto hay", espa, "espacios")
## En el texto hay 21487 espacios
tabla_porcen <- data.frame(letras=letters, porcentaje=porcen )
tabla_porcen
##    letras  porcentaje
## 1       a  7.64242169
## 2       b  1.50199706
## 3       c  1.94134959
## 4       d  3.99621610
## 5       e 12.41013244
## 6       f  1.90666386
## 7       g  1.75425688
## 8       h  6.92243010
## 9       i  5.97960900
## 10      j  0.07778011
## 11      k  0.86609207
## 12      l  4.63107000
## 13      m  3.01240277
## 14      n  6.33487492
## 15      o  8.31826782
## 16      p  1.31175110
## 17      q  0.06621820
## 18      r  6.02795880
## 19      s  6.55455119
## 20      t  9.15598066
## 21      u  3.49064536
## 22      v  1.05423586
## 23      w  2.28715577
## 24      x  0.13138533
## 25      y  2.59302081
## 26      z  0.03153248

Teorema central del límite

Es bien conocido el teorema que abordamos en este ejercicio y solo para recordar, si \(X_1, X_2, \ldots\) es una sucesión de variables aleatorias independientes e identicamnete distribuidas, con media \(\mu\) y varianza finita \(\sigma^2\), la función de distribución de la variable aleatoria \(Z\) descrita por

\[Z = \frac{(X_1 + \ldots + X_n) - n\mu}{\sqrt{n\sigma^{2}}}\]

tiende a la función de distribución normal estándar cuando \(n \to \infty\). Entonces tu objetivo será crear una función que, de acuerdo a una distribución (pueden ser tantas como conozcas) que sean adaptables a las condiciones del teorema, incluyendo como parámetros de la función media y la varianza de dicha distribución y una \(n\), se creen simulaciones de dicha distribución, al igual que la v.a. \(Z\). Finalmente se tiene que dar las gráficas de probabilidad acumulada y de densidad correspondiente (puedes guardar todo en una lista). Al final, con un \(n\) grande dado en la función se debería ver una “comprobación visual de dicho teorema”.

Para empezar vamos a ver la distribución poisson.

Para el caso en el que X se distribuya Poisson con parámetro \(\lambda\), tenemos que ver que:

  1. \(\mathbb{E}(X) = \lambda\)
  2. \(\mathbb{Var}(X) = \lambda\)

Vamos a simular para \(n=300\) y \(\lambda = 2\)

En este caos \(\mathbb{E}(X) = 2\) y \(\mathbb{Var}(X) = 2\) Para generar la simulación utilizamos y modoficamos un codigo creado por el profesor Luis Rincón

k <- 25000 #tamaño de las simulaciones 
v <- rep(0,k)
n <- 300 # el numero de sumandos

lam <- 2 #parámetro de la poisson
for(i in 1:k){
  x <- rpois(n,lam)
  v[i] <- sum(x)
}

Realizamos ahora la estandarización

v <- (v-n*lam)/sqrt(n*lam)

Para darnos una idea veamos la funcion de densidad y probabilidad acumulada de una Poisson

plot(dpois(0:10,2),type="b",xlab="x",ylab="P(X=x)",main="Función de densida Poisson(2)", pch=19, col="black")

plot(ppois(0:10,2),type="b",xlab="x",ylab="P(X=x)",main="Probabilidad acumulada Poisson(2)", pch=19, col="black")

Ahora hagamos el histograma de la estandarización para visualizar si efectivamente tinene la forma de una normal estandar

hist(v, main = "Límite central", xlim = c(-5,5))

Vemos que si se parce a una una distribución normal con \(\mu = 0\) y \(\sigma^2 =1\), es decir, es una normal estandar.