Paradoja del cumpleaƱos

Se desea determinar la probabilidad de que dos personas en un salĆ³n cumplan el mismo dĆ­a. Se tienen \(n\) personas, los aƱos bisiestos no son contados ni se admiten las personas gemelas; ademĆ”s de que los posibles 365 cumpleaƱos tienen la misma probabilidad de ocurrir.

Se tienen las siguientes 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 tĆŗ. \[P_1=\mathbb{P}[\text{2 personas cumplen aƱos el mismo dĆ­a}]=\begin{cases} 1āˆ’\frac{365!}{365^n(365āˆ’n)!}~~~~~~1ā‰¤nā‰¤365\\1~~~~~~~~~~~~~~~~~~~~~~~~~~~~~n>365 \end{cases}\]

\[P_2=\mathbb{P}[\text{Otra persona cumple aƱos el mismo dia que tĆŗ}]=1āˆ’\left(\frac{364}{365}\right)^n\]

1. Crea una funciĆ³n que, de acuerdo a una \(n\) vĆ”lida, determine ambas probabilidades.

#funcion que calcula ambas probabilidades, que dos personas cumplan el mismo dia y 
#que otra persona cumpla el mismo dia que tĆŗ
#recibe como parametro la cantidad de personas
proba_cumple<-function(n){
  p1=0
  if(n>365){
    p1=1
  }else if(n>=1 & n<=365){
    p1=1-prod(seq(365-n+1,365))/(365^n)
  }else{
    print("El nĆŗmero ingresado no es vĆ”lido")
  }

  p2=1-(364/365)^n
  return(list(p1=p1,p2=p2))
}

2. Crea una grĆ”fica donde se tengan la distribuciĆ³n de cada una de las probabilidades y determina 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 personas cumpla el mismo dĆ­a que tĆŗ.

library(ggplot2)

#La funcion recibe la cantidad de personas
graficar_cumple<-function(n){
  personas<-seq(1,n)
  probabilidad1<-seq(1,n)
  probabilidad2<-seq(1,n)
  for(i in 1:n){
    p=proba_cumple(i)
    probabilidad1[i]=p$p1
    probabilidad2[i]=p$p2
  }
  
  datos<-data.frame(personas, probabilidad1, probabilidad2)
  
  ggplot(datos) + 
    geom_line(aes(x=personas, y=probabilidad1), colour="red") + 
    geom_line(aes(x=personas, y=probabilidad2), colour="blue") +
    labs(title="Paradoja del cumpleaƱos", 
         y="Probabilidad", 
         x="Cantidad de personas")+
    theme(plot.title = element_text(hjust = 0.5))+
    annotate("text", x=n+n/10, y=0.9*probabilidad1[n], label= "P1", colour = "red")+
    annotate("text", x=n+n/10, y=0.5*probabilidad2[n], label= "P2", colour = "blue")
    
}

graficar_cumple(2000)

Se puede observar en la grĆ”fica que las funciones intersectan en algĆŗn valor pequeƱo de \(n\),y ademĆ”s parecen aproximarse entre sĆ­ cuando \(n\to \infty\), reducimos la cantidad de personas, sea \(n=5\):

Se puede observar que las probabilidades aparentemente coinciden para \(n=3\), calculando ambas probabilidades tenemos que:

## [1] "La probabilidad de que en un grupo de 3 personas(tĆŗ incluido), 2 cumplan aƱos el mismo dĆ­a es 0.008 y la probabilidad de que una de ellas cumpla aƱos el mismo dĆ­a que tĆŗ es 0.008"

TambiĆ©n para algĆŗn valor grande de \(n\) estas probabilidades se asemejan, por ejemplo para \(n=2000\)

## [1] "La probabilidad de que en un grupo de 2000 personas(tĆŗ incluido), 2 cumplan aƱos el mismo dĆ­a es 1 y la probabilidad de que una de ellas cumpla aƱos el mismo dĆ­a que tĆŗ es 0.996"

Entonces en un grupo de 3 personas, y para grupos muy numerosos(de 2000 personas o mƔs) las probabilidades \(P1\) y \(P2\) son iguales.

RelaciĆ³n Fibonacci-Eigen (vectores/valores)

Los nĆŗmeros de Fibonacci quedan representados por la ecuaciĆ³n recursiva \(Fn=Fnāˆ’1+Fnāˆ’2\) y de una manera muy sencilla se puede 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}\]

Crea una funciĆ³n para obtener el n-Ć©simo nĆŗmero de Fibonacci.

#funcion que calcula el n-esimo numero de Fibonacci
fibonacci<-function(n){
  #se definen los 2 primeros tƩrminos de la sucecion
  f0=0
  f1=1
  if(n==0){
    f=f0
  }else if(n==1){
    f=f1
    #se definen los siguientes termimos de forma iterativa
  }else if(n>=2){
    j=n-1
    for(i in 1:j){
      f=f1+f0
      f0=f1
      f1=f
    }  
  }
  return(f) 
}

1.Determina mediante el uso de R el eigen valor positivo correspondiente a dicha matriz (es decir, el famoso nĆŗmero Ć”ureo o nĆŗmero de oro).

#buscamos los eigen valores de la matriz y seleccionamos el positivo  
m<-matrix(c(1,1,1,0), nrow=2)

#los eigen valores ya vienen ordenados de manera descendente asi que el eigenvalor
#positivo es el primero
valor<-eigen(m)$values[1]
valor
## [1] 1.618034

Entonces para la matriz \(\begin{pmatrix} 1 & 1\\ 1 & 0 \end{pmatrix}\) el eigen valor positivo es \(1.618034\) el cuĆ”l coincide con \(\varphi=\frac{1+\sqrt{5}}{2}\) llamado nĆŗmero Ć”ureo.

2.Crea una grĆ”fica, para un \(n\) que desees, donde cada punto corresponda a \((F_n,F_{nāˆ’1})\). Dichos puntos deben ser de color negro.

#graficamos los puntos (Fn,Fn-1), usando la siguiente funcion:
graficar_fib<-function(n){
  sec1<-vector()
  sec2<-vector()
  #llenamos los vectores vacios con las parejas Fn,Fn-1
  for(i in 1:n){
    sec1[i]=fibonacci(i)
    sec2[i]=fibonacci(i-1)
  }
  
  datos<-data.frame(sec1, sec2)
  
  p<-ggplot(datos, aes(x=sec1, y=sec2)) + 
    geom_point() + 
    labs(title="Sucesion de Fibonacci", y="Fn-1", x="Fn")+
    theme(plot.title = element_text(hjust = 0.5))
  return(p)
}

n=10
graficar_fib(n)

3.En la misma grƔfica coloca la recta sobre la que pasa el eigen vector correspondiente al eigen valor del punto 1.

#graficamos la recta que tiene como vector director, el eigen vector asociado al
#eigen valor positivo
n=10

vec<-eigen(m)$vectors[,1]
h<-vector()
v<-vector()
for(i in 1:n){
  h[i]=(i-1)*(fibonacci(n)/(n-1))
  v[i]=h[i]*(vec[2]/vec[1])
}

graficar_fib(n)+geom_line(aes(x=h, y=v), colour="blue")

4.Elige algĆŗn punto de los graficados en el punto dos y multiplĆ­calo por el eigen valor del punto 1 y grafĆ­calo en color rojo ĀæQuĆ© sucediĆ³?

Tomamos el punto \((F_8,F_7)=(21,13)\)

punto<-c(fibonacci(n-2)*valor,fibonacci(n-2-1)*valor)
#graficamos el punto
graficar_fib(n)+geom_line(aes(x=h, y=v), colour="blue")+
  geom_point(aes(punto[1], punto[2]), colour= "red")

Resulta que al multiplicarlo por el eigen valor \(\varphi\), se obtiene que \(\varphi\cdot(F_8,F_7)=(\varphi\cdot F_8,\varphi\cdot F_7)\approx(F_9,F_8)\)

5.ĀæQuĆ© concluyes de todo esto?

A medida que \(n\) crece la razĆ³n entre los nĆŗmeros \(F_n\) y \(F_{n-1}\) se va aproximando a \(\varphi\), es debido a esto que los puntos \((F_n,F_{n-1})\) a medida que n crece parecen ajustarse mĆ”s a la recta cuyo vector director es el eigen vector asociado al eigen valor \(\varphi\), en otras palabras:

\[\displaystyle\lim_{n \to \infty}{\frac{F_n}{F_{n-1}}}=\varphi\]

Iteraciones

1. Derivada

Para comprobar determina si la derivada de \(2x^2\) en algĆŗn punto se aproxima con tu funciĆ³n. Para aproximar la derivada usamos la expresion: \[f'(x_o)=\frac{f(x_o+h)-f(x_0)}{h}\]

library(latex2exp)
library(ggplot2)
library(rSymPy)
library(knitr)
#funcion que calcula la derivada exacta, recibe a f como una cadena y el valor en
#el que serĆ” evaluada la derivada
derivada <- function(f,punto){
  x<-Var('x')
  caract <- paste("expr =", f) 
  eq <- sympy(caract)
  deriv <- sympy("d_expr = diff(expr, x, 1)")
  ev <- sympy(paste("d_expr.subs(x,", toString(punto),")"))
  return(ev)
}
#derivada aproximada
deriva_aprox<-function(punto,h){
 f = function(x) 2*x^2 #aqui se cambia la expresion dependiendo de la funcion 
 d=(f(punto+h)-f(punto))/h
 return(d)
}

#funcion que grafica la derivada exacta y la aproximada
graf<-function(f){
  x <- seq(-5, 5, length=100)
  y<-vector()
  y2<-vector()
  for(i in 1:length(x)){
    y[i]=as.numeric(derivada(f,x[i]))
    y2[i]=deriva_aprox(x[i],0.1)
  }
  p1<-ggplot(data.frame(x,y))+
    geom_line(aes(x,y),color="steelblue", size=1, linetype="dashed")+
    geom_line(aes(x,y2),color="red", size=1, linetype="dotted")+
    ggtitle("Derivada") +
    labs(y= "f'(x)", x = "x") +
    theme(plot.title = element_text(hjust = 0.5)) +
    annotate("text", x=4, y=4, label= "derivada", colour = "steelblue")+
    annotate("text", x=4, y=3, label= "aproximaciĆ³n", colour = "red") 
  return(p1)
}

Evaluamos la funcion \(2x^2\) y graficamos su derivada exacta y aproximada:

graf("2*x**2")

Podemos ver que con una \(h=0.1\) nuestra derivada aproximada, se asemeja a la derivada real.

2. Integral

Puedes usar funciones positivas para comprobar tu funciĆ³n utilizando la interpretaciĆ³n de la integral. Como tal no existe (o a estas alturas no conocemos) un mĆ©todo para resolver integrales de forma iterativa, asĆ­ que generaremos una funciĆ³n que permita aproximar una integral usando integraciĆ³n numĆ©rica. Usamos el mĆ©todo del trapecio: \[\int_a^bf(x)dx\approx\frac{h}{2}\left(f(a)+f(b)+2\sum_{i=1}^{n-1}f(x_{i})\right) \\ \text{con}~~~h=\frac{b-a}{n},~~~x_i=a+ih,~~~i=1,2,3,...,n\]

#funcion que calcula la integral exacta, recibe la funcion como parametro, los limites de
#integracion y devuelve la integral evaluada
integral<-function(f,a,b){
  return(integrate(f,a,b))
}

#integral aproximada
#usamos integracion numerica, recibe como parametro la funcion
# los limites de integracion y la cantidad de divisiones que se
# haran en el intervalo.
#n debe ser par
integral_aprox<-function(f,a,b,n){
  h=(b-a)/n
  i=seq(1,n-1)
  inte<-h/2*(f(a)+f(b)+2*sum(f(a+i*h)))
  return(inte)
}

Comparamos los resultados para la siguiente integral: \[\int_0^2(x^3-4x^2+4x+1)~dx\]

La integral exacta es:

f=function(x) x**3-4*x**2+4*x+1
integral(f,0,2)
## 3.333333 with absolute error < 3.7e-14

Ahora, variando la cantidad de subintervalos(\(n\)) aproximemos la integral con la regla de Simpson.

library(knitr)

f=function(x) x**3-4*x**2+4*x+1
m<-seq(4,100,5)

aproximacion<-vector()
for(i in 1:length(m)){
  aproximacion[i]=integral_aprox(f,0,2,m[i])
}

library(knitr) 
kable(data.frame(n = m,Aproximacion=aproximacion))
n Aproximacion
4 3.250000
9 3.316872
14 3.326531
19 3.329640
24 3.331018
29 3.331748
34 3.332180
39 3.332457
44 3.332645
49 3.332778
54 3.332876
59 3.332950
64 3.333008
69 3.333053
74 3.333090
79 3.333120
84 3.333144
89 3.333165
94 3.333182
99 3.333197

Se puede observar que para valores grandes de n, la aproximaciĆ³n de la integral se asemeja mĆ”s a la integral real.

3. PerĆ­metro de una circunferencia

Para un polĆ­gono regular de \(n\) lados inscrito en una circunferencia de radio \(r\) se sabe que la medida de uno de sus lados estĆ” dada por: \[L=2r\cdot\sin\left( \frac{\pi}{n}\right)\]

Y su perimetro es: \[P=n\cdot L\]

Tambien recordemos que dada la medida del lado de un poligono regular de \(n\) lados inscrito en una circunferencia de radio \(r\), sea \(L_n\) esa medida, la medida del lado de un poligono de \(2n\) lados inscrito en la misma circunferencia es: \[L_{2n}=\sqrt{2r^2-r\sqrt{4r^2-L_{n}^2}}\]

Crea una funciĆ³n, que de acuerdo a un numero de lados n se vaya acercando al perĆ­metro de una circunferencia. Al final puedes comprobar tus resultados con la formula ya conocida.

Necesitamos una funciĆ³n que vaya mejorando el resultado con cada iteraciĆ³n y que cada paso dependa del paso anterior.

#la funcion recibe el radio, numero de lados y cantidad de iteraciones
perimetro<-function(r,n,it){
  l=2*r*sin(pi/n)
  a<-vector()
  b<-vector()
  a[1]=l*n
  b[1]=n
  for(j in 1:it){
    #usamos un metodo iterativo para calcular el valor del lado de un poligono de 2n lados
    #a partir de un poligono de n lados
    l=sqrt(2*r^2-r*sqrt(4*r^2-l^2))
    n=2*n
    a[j+1]=l*n
    b[j+1]=n
  }
  return(list(a=a,b=b))
}

Comprobamos que la funciĆ³n sirve usando la formula usual para el perimetro de la circunferencia de radio \(r\): \[P=2\cdot\pi\cdot r\]

Para una circunferencia de radio \(r=6\), su perimetro es: \[P=2\cdot\pi\cdot 6\approx37.69911\]

Usando la funciĆ³n, empezando con un poligono de 4 lados \(n=4\), realizando 15 iteraciones, se tiene que:

library(knitr) 
p<-perimetro(6,4,15)
kable(data.frame(n = p$b,Perimetro = p$a))
n Perimetro
4 33.94113
8 36.73761
16 37.45734
32 37.63858
64 37.68397
128 37.69533
256 37.69817
512 37.69888
1024 37.69905
2048 37.69910
4096 37.69911
8192 37.69911
16384 37.69911
32768 37.69911
65536 37.69911
131072 37.69911

Se observa que mientras n crece, el perĆ­metro del polĆ­gono se aproxima bastante al de la circunferencia.

Procesamiento de textos

De acuerdo a un archivo .txt determinar lo siguiente:
Cantidad de letras.
Cantidad de vocales.
Cantidad de espacios.
Porcentaje que representa cada letra en el texto.

library(tidyverse)
#Coloca aqui la direccion o nombre del texto que desea usar
direccion<-"Texto.txt"
texto<-read.table (direccion, sep="\r",encoding="UTF-8")
texto = toupper(texto$V1)

#transformamos cada renglon del texto en una lista de letras
texto_separado = strsplit(texto, split="")

#Deshacemos esa lista y tenemos el data.frame
texto_col = as.character(unlist(texto_separado))
texto_col = data.frame(texto_col)
names(texto_col) = c("V1")

#obtenemos la cantidad de espacios antes de eliminar caracteres innecesarios
df<-data.frame(table(texto_col))
for(i in 1:length(df$texto_col)){
  if(df$texto_col[i]==" "){
    espacios=df$Freq[i]
  }
}

#Eliminamos signos de puntuacion, caracteres y numeros
texto_col$V1 = sub("([[:space:]])","",texto_col$V1)
texto_col$V1 = sub("([[:digit:]])","",texto_col$V1)
texto_col$V1 = sub("([[:punct:]])","",texto_col$V1)

#quitamos los acentos de las vocales
for(i in 1:length(texto_col[[1]])){
  if(texto_col[[1]][i]=="Ɓ"){
    texto_col[[1]][i]="A"
  }else if(texto_col[[1]][i]=="Ɖ"){
    texto_col[[1]][i]="E"
  }else if(texto_col[[1]][i]=="ƍ"){
    texto_col[[1]][i]="I"
  }else if(texto_col[[1]][i]=="Ɠ"){
    texto_col[[1]][i]="O"
  }else if(texto_col[[1]][i]=="ƚ"){
    texto_col[[1]][i]="U"
  }
}
  
#contamos cuantas letras y vocales hay
tabla=rev(table(texto_col))
letras=sum(tabla[1:26])
vocales=tabla[5]+tabla[11]+tabla[18]+tabla[22]+tabla[26]

#obtenemos el porcentaje que representa cada letra en el texto
frecuencias=round(100*tabla[1:26]/letras, digits=2)

t<-data.frame(frecuencias)
t$texto_col <- unlist(t$texto_col)
t$Freq <- unlist(t$Freq)

library(ggplot2)

letra<-t$texto_col
porcentaje<-t$Freq
#graficamos
p<-ggplot(t, aes(letra, porcentaje)) +
  geom_bar(stat = "identity", color = "black", fill = "#87CEFA") +
  geom_text(aes(hjust = 1.3, label = porcentaje)) + 
  coord_flip() + 
  labs(title = "Porcentaje de cada letra en el texto",  x = "Letra", y = "Porcentaje")+
  theme(plot.title = element_text(hjust = 0.5))
library(knitr) 
m<-c("Espacios","Letras","vocales")
o=c(espacios,letras,vocales)
kable(data.frame(Objetos = m,Cantidad=o))
Objetos Cantidad
Espacios 612
Letras 3084
vocales 1417 
p

Teorema central del lĆ­mite

Para \(X_1,X_2,X_3ā€¦\) una sucesiĆ³n de variables aleatorias independientes e idĆ©nticamente 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+ā‹Æ+X_n)āˆ’n\mu}{\sqrt{n\sigma^2}}=\frac{\bar X-\mu}{\frac{\sigma}{\sqrt{n}}}\] tiende a la funciĆ³n de distribuciĆ³n normal estĆ”ndar cuando \(n \to \infty\).

Crear una funciĆ³n que, de acuerdo a una distribuciĆ³n que sea adaptable a las condiciones del teorema, incluyendo como parĆ”metros de la funciĆ³n la 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.

#media: media de la muestra
#varianza: varianza de la muestra
#n: tamaƱo de la muestra
#m: numero de simulaciones
TCL<-function(media, varianza, n, m){
  z<-vector()
  for(i in 1:m){
    #aqui se cambia rchisq por la distribucion que vayas a usar y sus respectivos parametros
    muestra<-rchisq(n, lamda) 
    z[i]<-(mean(muestra)-media)/sqrt(varianza/n)
  }
  z=data.frame(z)
  
  #grĆ”fica de densidad de dicha distribuciĆ³n
  x <- seq(0, 10, length=10)
  y <- dchisq(x, lamda)# esta se cambia dependiendo la distribucion y sus parametros
  p1<-ggplot(data.frame(x,y),aes(x,y))+
    geom_line()+
    ggtitle("FunciĆ³n de Densidad") +
    labs(y= "f(x)", x = "x") +
    theme(plot.title = element_text(hjust = 0.5))
  
  #grĆ”fica de prob. acum de dicha distribuciĆ³n
  y <- pchisq(x, lamda)# esta se cambia dependiendo la distribucion y sus parametros
  p2<-ggplot(data.frame(x,y),aes(x,y))+
    geom_line()+
    ggtitle("FunciĆ³n de DistribuciĆ³n") +
    labs(y= "F(x)", x = "x") +
    theme(plot.title = element_text(hjust = 0.5))  
  p<-grid.arrange(p1,p2, ncol=2, nrow=1)
  
  g<-ggplot(z, aes(x = z)) + 
    geom_histogram(fill = "blue",alpha = .20, binwidth=0.2, colour = "black", aes(y = ..density..)) + 
    geom_density(colour="blue", size=1) +
    stat_function(fun = dnorm, colour = "red", size= 1, args = list(mean = 0, sd = 1))+
    xlab("x") + 
    ylab("Densidad")+
    annotate("text", x=3, y=0.4, label= TeX('$N(0,1)$'), colour = "red")+
    annotate("text", x=3, y=0.35, label= TeX('$f(Z)$'), colour = "blue")
  return(list(p=p,g=g))
  
}

Finalmente se tiene que dar las grĆ”ficas de probabilidad acumulada y de densidad correspondiente . Al final, con un n grande dado en la funciĆ³n, se deberĆ­a poder ver una ā€œcomprobaciĆ³n visual de dicho teoremaā€.

Para ejemplificar, la muestra \(X_1, X_2,....\) tendrĆ” una distribuciĆ³n \(\chi^2(\lambda)\) con \(\lambda=3\), para esta distribuciĆ³n se tiene que: \[\mu=\lambda~~~~~~y~~~~~~~\sigma^2=2\lambda\] A continuaciĆ³n se muestran las grĆ”ficas y la simulaciĆ³n:

library(latex2exp)
library(ggplot2)
require(gridExtra)

#tamaƱo de la muestra
n=100
#definimos a que distribucion pertenecerĆ” la muestra, en este caso pertenece a la distribucion
#Xi cuadrada con parametro lambda=3
lamda=3
#para esta distribucion en particular la media y la varianza son:
media=lamda
varianza=2*lamda
TCL(media, varianza, n, 10000)$g

Se puede observar que para valores grandes de n, la distrbuciĆ³n de \(Z\)(azul) tiende a la Normal Estandar(rojo).