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}\]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.
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}\]
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.
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:
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
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:
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.