SIMULACIÓN ESTOCÁSTICA
TAREA 1
Alumna: Andrea Elizabeth Arredondo Sánchez
1. Verificar que el generador estándar mínimo satisface la proposicion 1
Gen.estandar.min <- function(n,a,m,c, sem=as.numeric(Sys.time())){
U=numeric(n)
U[1] <- sem
for(i in 1:(n-1)){
U[i+1]<- (a*U[i]+c)%%m
}
return(U/m) #genero una secuencia de numeros pseudoaleatorios (aproximac. de uniformes)
}
a=7^5; m= 2^31-1; c=0
p1<- Gen.estandar.min(200,a,m,c) #tomo una muestra de tamaño n=200
summary(p1) #verifico que sean Uniformes(0,1)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.001864 0.220450 0.479852 0.486603 0.730848 0.998973Al implementar el algoritmo ‘Gen.estandar.min’ en la consola se obtuvo una muestra de 200 observaciones de una variable aleatoria aproximada a la Uniforme en (0,1)
Luego traté de corroborar que el generador mínimo estandar efectivamente satisface la Proposición 2.1 de las notas; sin embargo como c=0 eso genera que la primer condición no se cumpla. Es decir, el máximo común divisor mcd(c,m)=m= 2^31-1 y es distinto de la unidad.
2. Hacer pruebas estadísticas a los numeros pseudoaleatorios generados en el ejercicio anterior
rachas <- factor()
datos_dicotomicos <- function(datos){ #genero una funcion que convierta datos de U(0,1) en dicotomicos
for (i in 1:length(datos)){
a<-0
if(datos[i]>0.5) {a<-1} #defino punto de corte
rachas<-c(rachas,a)
}
return(rachas)
}
(p2<-datos_dicotomicos(p1)) #asigno rachas a valores generados por p1## [1] 1 1 0 0 0 1 0 1 0 1 1 0 1 0 0 0 1 0 1 1 0 1 1 0 1 1 0 0 0 0 1 1 0 0 1
## [36] 0 1 0 0 1 0 1 1 1 0 0 1 1 0 1 0 1 0 1 0 1 1 1 0 0 1 0 0 1 1 0 0 0 1 1
## [71] 1 1 1 0 0 1 1 1 0 0 1 0 1 0 0 1 0 1 1 0 0 1 0 1 0 1 1 0 1 0 0 0 0 1 1
## [106] 1 0 1 0 1 0 0 1 1 1 0 0 1 1 1 1 0 0 1 0 0 1 1 0 0 1 0 0 1 1 1 1 0 1 1
## [141] 1 1 1 1 0 0 0 1 1 1 0 1 1 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 1
## [176] 1 1 1 0 0 0 1 0 0 0 1 1 0 0 1 1 1 0 0 0 0 1 1 0 0##Aplico prueba de rachas para ver si los datos conservan su aleatoriedad
library(adehabitatLT)
wawotest(p2) ## a ea va za p
## -8.0428171 -1.0000000 197.9983913 -0.5005133 0.6916431## Aplico prueba de Anderson-Darling para ver que distribuyen U(0,1)
library(ADGofTest)
ad.test(p2,punif,0,1) ##
## Anderson-Darling GoF Test
##
## data: p2 and punif
## AD = Inf, p-value = 3e-06
## alternative hypothesis: NAPrimero, genere la función ‘datos_dicotomicos’ que convierte los datos que se le ingresen en dicotómicos para poder aplicar las pruebas estadísticas
Como resultado de la prueba de rachas obtuve un p-value = 0.07147075 > 0.5 por lo cual acepto H0 que indica que los datos si son aleatorios y de la prueba de A-D obtuve un p-value = 3e-06 por lo cual también se corrobora que los datos distibuyen U(0,1)
3. Comparar la funcion de distribucion empírica con la funcion de distribucion teorica de U(0,1) de manera gráfica
dist.empirica <- ecdf(p1)
x=rnorm(200) #muestra de tamaño n=200
qqplot=qqnorm(x)
plot(dist.empirica,col="blue",lwd=3)
points(qqplot$y,qqplot$x,col="red")
Al observar la gráfica se verifica que son muy parecidas ambas curvas. Además,dichas curvas son muy parecidas a la identidad y eso quiere decir que x% de los valores de la muestra son menores o iguales que x
4. GENERADOR de Wichmann y Hill (1982)
x<-c();y<-c();z<-c()
wichman_hill <- function(n,a,b,c){
if(a>0 & a<30000 & b>0 & b<30000 & c>0 & c<30000){ #defino intervalo de las variables
v <- c(30269,30307,30323)
u <- rep(0,n)
u[1]<-(sum(c(a,b,c)/v)) %% 1
x[1]=a; y[1]=b; z[1]=c #los valores iniciales (semillas) otorgados por el usuario
for(i in 1:(n-1)){
if(x[i]>0){x[i+1] <- (171*x[i])%%177 - (2*x[i]/177)}
else{x[i+1]=x[i]+30269}
if(y[i]>0){y[i+1] <- (172*y[i])%%176 - (35*y[i]/176)}
else{y[i+1]=y[i]+30307}
if(z[i]>0){z[i+1] <- (170*z[i])%%178 - (63*x[i]/178)}
else{z[i+1]=z[i]+30323}
u[i+1] <- ((x[i+1]/30269)+(y[i+1]/30307)+(z[i+1]/30323))%%1
}
} else{print("Los números recibidos por la función no están en el intervalo (0,30000)")}
return(u) #regresa los numeros generados
}
p4<-wichman_hill(500,0.5,7000,18065)
evaluar<-datos_dicotomicos(p4) #hago dicotomicos para hacer pruebas estadisticas
ad.test(evaluar,punif,0,1) #p-value=1.2e-06 entonces si distribuyen U(0,1)##
## Anderson-Darling GoF Test
##
## data: evaluar and punif
## AD = Inf, p-value = 1.2e-06
## alternative hypothesis: NAlibrary(randtests)
runs.test(evaluar,threshold=0.5) #rechazo aleatoriedad (p-value < 2.2e-16)##
## Runs Test
##
## data: evaluar
## statistic = -22.22, runs = 2, n1 = 53, n2 = 447, n = 500, p-value
## < 2.2e-16
## alternative hypothesis: nonrandomnessImplementé el algoritmo propuesto por wichman & Hill para generar números pseudoaleatorios realizando 500 iteraciones dados los valores iniciales a=0.5, b=7000 y c=18065 que están dentro del intervalo (0,30000) y luego apliqué el test de rachas y el de bondad de ajuste.
Los resultados que obtuve es que los datos si distribuyen U(0,1) y rechazo aleatoriedad con punto de corte en p=0.5 y eso se ve claramente por la manera en la que fueron arrojados los valores luego de hacerlos dicotómicos (en la variable evaluar) por lo cual el generador no es confiable
5. Considere una variable aleatoria X con funcion de densidad f(x) = 2x si xE(0,1). Escribir e implementar un algoritmo para generar valores de dicha variable.
generador_continua<- function(n){
V<-numeric(n)
for(i in 1:n){
U<-runif(1)
V[i]<- (U)^(1/2)}
return(V)
}
p5<-generador_continua(50)
summary(p5)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.1340 0.4471 0.7385 0.6559 0.8328 0.9870Por ser una función continua decidí usar la función inversa para obtener 50 valores aleatorios en (0,1) y debido al ‘summary’ mostrado anteriormente se puede corroborar que los datos permanecen en dicho intervalo.
6. Considere una variable aleatoria X, con soporte en el conjunto {1,2,3,4,5} y con distribucion {1/3; 1/5; 1/6; 1/6; 4/30}. Implementar un algoritmo para generar una muestra de 1, 000 elementos de dicha distribucion
Tinv_discretas <- function(n,D){
if(sum(D)!=1)stop("El vector de densidades debe sumar 1")
muestra<-numeric(n)
for(i in 1:n){
k<-0
U<-runif(1) #genero 1000 uniformes(0,1)
F<-0
while(F<U){
F<-F +D[k+1]
if(F<U){
k<-k+1
}else{F<-1}
}
muestra[i]<-k
}
return(muestra)
}
p6<-Tinv_discretas(1000,c(1/3,1/5,1/6,1/6,4/30))
table(p6)## p6
## 0 1 2 3 4
## 351 198 151 163 137Con el uso del algoritmo ‘Tinv_discretas’ de la Transformada Inversa para el caso de las variables discretas generé una muestra de tamaño n=1000 según la distribución que se indicó y obtuve los resultados anteriores.
7. Implementar algoritmos para generar numeros provenientes de las siguientes variables aleatorias:
a) Uniforme discreta con soporte en los primeros 10 naturales
p7_a<- Tinv_discretas(1000,c(rep(1/10,10)))
table(p7_a)## p7_a
## 0 1 2 3 4 5 6 7 8 9
## 100 86 84 111 100 94 113 100 116 96Usando el mismo algoritmo que en el ejercicio anterior obtengo los resultados enunciados para cada valor que toma la Uniforme
b)Binomial con n = 10 y p = 1/2
G.bin.inv.recursiva<-function(n,p){ #generador de observaciones Bin(n,p) por inversa generalizada
V<-numeric(n)
for(j in 1:n){
c<- (p/(1-p))
pr<- (1-p)^n
F<- pr
i<- 0
U<- runif(1)
if(F<U){
while(F<U){
i<-i+1
pr<- c*((n-i+1)/i)*pr
F<-F+pr
}
}
V[j]<-i
}
return(V)
}
(p7_b<- G.bin.inv.recursiva(40,1/2))## [1] 23 15 19 19 21 17 20 19 17 21 21 21 21 22 22 20 26 20 24 23 19 20 24
## [24] 24 18 16 16 23 23 22 20 19 28 28 23 22 16 17 24 18Con la implementación de las probabilidades iterativas de la variable Binomial se pudieron simular n=40 valores de X con una probabilidad de p=1/2 con el algoritmo ‘G.bin.inv.recursiva’ y los valores obtenidos son los que se guardaron en la variable p7_b
c) Poisson(9)
G.poi.inv.recursiva<-function(n,lambda){
V<-numeric(n)
for(j in 1:n){
p<-exp(-lambda)
F<-p
i<-0
U<-runif(1)
if(F<U){
while(F<U){
i<-i+1
p<-p*lambda/i
F<-F+p
}
}
V[j]<-i
}
return(V)
}
p7_c<- G.poi.inv.recursiva(100,9) ; p7_c ## [1] 8 8 11 15 8 8 10 11 8 10 10 9 13 5 9 16 5 6 13 12 6 9 10
## [24] 7 16 12 7 12 8 9 12 11 9 6 13 9 10 10 14 13 9 14 8 8 16 6
## [47] 11 7 10 9 9 12 15 11 4 4 8 13 10 6 9 12 8 11 5 6 10 9 9
## [70] 11 11 12 9 4 11 6 9 5 6 13 9 8 7 9 9 6 7 6 5 6 5 6
## [93] 12 10 12 6 6 9 5 8plot(rpois(100,9));plot(p7_c)
Usando el generador de observaciones Poisson por medio de la inversa generalizada con las probabilidades iterativas de dicha distribución simule una muestra de tamaño n=100 con lambda=9 y luego comparé graficamente dichos valores con la función teórica de probabilidad y en efecto se parecen mucho los valores simulados a los datos reales
d)Geometrica con un parametro adecuado para tener un numero esperado de iteraciones menor a 3
G.geom.inv <-function(n,p){
V<-numeric(n)
for(j in 1:n){
for(i in 1:400){
U<-runif(1)
if(p<U){i<-i+1}
else{
break
}
}
V[j]<-i
}
return(V)
}
p7_d<- G.geom.inv(200,1/4)
table(p7_d)## p7_d
## 1 2 3 4 5 6 7 8 9 10 11 12 14 16
## 60 34 26 18 16 10 12 4 8 4 1 4 2 1Dado que X ~ Geometrica (1/4) entonces el número esperado de iteraciones resulta menor a 3 y dado que X = número esperado de ensayos necesarios para que ocurra por PRIMERA VEZ un éxito entonces generé una función llamada ‘G.geom.inv’ cuya finalidad es simular n = 200 valores aleatorios que distribuyan Geom(1/4)
8. Describir un algoritmo para generar la distribucion Binom utilizando el algoritmo de la distribuion Bernoulli
Generador_binomial <- function(muestra,n,p){
X<-numeric(n)
r<-numeric(muestra)
for(j in 1:muestra) {
for(i in 1:n){ #realizamos tantas Blli como n
U<-runif(1)
if(U<=p){X[i]<-1}
else{X[i]<-0}
}
r[j]<-sum(X)
}
return(r)
}
p8<-Generador_binomial(40,50,1/2) #obtengo una muestra de 40 observaciones de una Bin(n,p)
table(p8)## p8
## 17 19 20 22 23 24 25 26 27 28 29 30 31
## 1 1 3 1 3 3 7 6 6 4 2 2 1Implementé la función ‘Generador_binomial(40,50,1/2)’ generando 40 valores que distribuyen Binom(50,1/2) y obtuve los resultados anteriores
9. Generar una muestra aleatoria de una distribucion Bin(1000, 1/2) utilizando la definicion de funcion inversa generalizada.
p9<- G.bin.inv.recursiva(1000,1/2)
summary(p9)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 445.0 490.0 501.0 500.5 511.0 558.0Implementando el algoritmo ‘G.bin.inv.recursiva’ utilizado y definido anteriormente obtengo una simulación de la variable Bin(1000,1/2)
10. Sea X ~ BinNeg(r; p)
a)Use la relacion de las distribuciones Binomial Negativa con la Geometrica para dar un algoritmo para generar valores de X e implementarlo.
Generador_binomial_negativa <- function(n,k,p){ #n>> grande
geom<- G.geom.inv(n,p)
if(length(geom)>k){ #cuento los éxitos
texto <- paste("la cantidad de ensayos de una Blli necesarios para que ocurra el k-esimo éxito (",k, ") son ",sum(geom[1:k]))
print(texto)
} else{
print("el tamaño de n debe de ser aún mayor para alcanzar el k-esimo éxito")
}
}
p10_a<-Generador_binomial_negativa(200,195,1/4)## [1] "la cantidad de ensayos de una Blli necesarios para que ocurra el k-esimo éxito ( 195 ) son 717"Con la implementación del algoritmo ‘G.geom.inv’ implementado en ejercicios anteriores generé la función ‘Generador_binomial_negativa’ para una variable BinNeg(k,p) que contará la cantidad de ensayos que tienen que llevarse a cabo para que ocurra el k-esimo éxito y el parámetro n que se solicita en la función ‘Generador_binomial_negativa’ indica la cantidad de veces que se va a generar una variable Geométrica. De preferencia n tiene que ser grande.
c)Usar la identidad recursiva anterior para dar un segundo algoritmo e implementarlo
G.bin.negativa.recursiva<-function(r,p){ #generador de observaciones Bin(n,p) por inversa generalizada
V<-numeric(r)
for(j in 1:r){
c<- (1-p)
pr<- ((1+p)^(-r))
F<- pr
i<- 0
U<- runif(1)
if(F<U){
i<-i+1
pr<- c*pr*((i-1)/(i-r))
F<-F+pr
}
V[j]<-i
}
return(V)
}
p10_c<- G.bin.negativa.recursiva(195,1/4)