Punto 1-

Solucion:

\(f(x)=20x(1-x)^3 .....0<x<1\) Funcion objetivo

\(g(x)=1 .....0<x<1\).funcion auxiliar

Algoritmo:

  1. Genere numeros aleatorios de la función auxiliar y tambien genere de una uniforme continua entre cero y uno \(Y \sim g\) ; \(U\sim U_{(0,1)}\)

  2. acepte \(x=y\), si \(u<=f(Y)/C*g(Y)\); $ C=máx{f(x)/g(x)}$

  3. En caso contrario, regrese al paso 1.

Recoredemos:
Si tenemos que \(g(x)=1\) \(0<x<1\) Funcion de densidad de probabilidad. Y si queremos generar numeros aleatorios de ella por el metodo de la inversion debemos conocer su \(G(x)\) que en este caso es X, y evaluar su inversa en una uniforme \(u{(0,1)}\). Por tanto para este caso, X=U donde U es una \(u{(0,1)}\) 
f=function(x){20*x*(1-x)^3} ##función objetivo
g=function(x){1} ###funcion auxiliar 

fg=function(x){(20*x*(1-x)^3)/1} #cociente entre las 2 funciones

Veamos la grafica de la funcion objetivo.

curve(f(x),from=0,to=1)  #curva de la funcion de densidad de beta desde o hasta 1

Ahora aplicamos el algoritmo:

x1<-0
Y<-runif(n=1,min = 0,max = 1) #genere de una uniforme#de la auxiliar, en este caso coincide
U<-runif(n=1,min = 0,max = 1)#genere de una uniforme
c1=optimize(fg,interval = c(0,1),maximum = TRUE)$objective ##Calculamos el C maximo
##criterio de decicion:

if(U<=(1/c1)*fg(U)){x1=Y}
x1
## [1] 0.1481245

Lo anterior solo esta corriedo una vez, Hagamoslo para muchas:

a. Para n=5.000

x<-0
nsim<-5000
k<-0 #es un contador
j<-0 # es un contador

X<- numeric(nsim) #donde guardara la información, o puede ir en vacio
c1=optimize(fg,interval = c(0,1),maximum = TRUE)$objective

while(k<nsim){ #para un control de que va a terminar en algun momento
  U<-runif(1)
  j<-j+1 #para ver cuantas veces ingresa en k<nsim
  Y<-runif(1) #de la funcion auxiliar g(x)
  #citerio de decision:
  
  if(U<=(1/c1)*fg(Y)){
    k<-k+1
    x[k]=Y
  }
}
#x por comodidad no los mostrare
j #cuantos entran al while
## [1] 10592
k
## [1] 5000
hist(x,freq = FALSE)
curve(f(x),from=0,to=1,add = TRUE)  #curva de la funcion de densidad de beta desde o hasta 1

Conclusion: Graficamente notamos que los datos simulados si probienen de la funcion f(x); los datos simulados se ajuastan a la curva teorica.

a. Para n=10.000

x<-0
nsim<-10000
k<-0 #es un contador
j<-0 # es un contador

X<- numeric(nsim) #donde guardara la información, o puede ir en vacio
c1=optimize(fg,interval = c(0,1),maximum = TRUE)$objective

while(k<nsim){ #para un control de que va a terminar en algun momento
  U<-runif(1)
  j<-j+1 #para ver cuantas veces ingresa en k<nsim
  Y<-runif(1) #de la funcion auxiliar g(x)
  #citerio de decision:
  
  if(U<=(1/c1)*fg(Y)){
    k<-k+1
    x[k]=Y
  }
}
#x por comodidad no los mostrare
j #cuantos entran al while
## [1] 21246
k
## [1] 10000
hist(x,freq = FALSE)
curve(f(x),from=0,to=1,add = TRUE)  #curva de la funcion de densidad de beta desde o hasta 1

Punto 2

Punto 3

Solución:

a. para n=5000

n <- 5000
p <- 0.4   # pongamos par esta probabilidad en particular 
u <- runif(n)
x <- as.integer(u >1-p) #(u > 0.6) is a logical vector
#x #estos son los x generadps pero no los mostrare en pantalla
plot(factor(x))

mean(x)  #para ver que la media estimada coincide con la teorica
## [1] 0.4028
var(x)
## [1] 0.2406003
ks.test(x,dbinom(x= 5000, size = 1, prob = 0.4)) ##prueba de bondad de aguste, como el p>0.05 lus datos si se ajustan a la distribucion bernulli
## Warning in ks.test(x, dbinom(x = 5000, size = 1, prob = 0.4)): cannot compute
## exact p-value with ties
## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  x and dbinom(x = 5000, size = 1, prob = 0.4)
## D = 0.4028, p-value = 0.9969
## alternative hypothesis: two-sided

b. para n=10000

n <- 10000
p <- 0.4   # pongamos par esta probabilidad en particular 
u <- runif(n)
x <- as.integer(u >1-p) #(u > 0.6) is a logical vector
#x
plot(factor(x))

mean(x)  #para ver que la media estimada coincide con la teorica
## [1] 0.3941
var(x)
## [1] 0.2388091
ks.test(x,dbinom(x= 10000, size = 1, prob = 0.4)) ##prueba de bondad de aguste, como el p>0.05 lus datos si se ajustan a la distribucion bernulli
## Warning in ks.test(x, dbinom(x = 10000, size = 1, prob = 0.4)): p-value will be
## approximate in the presence of ties
## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  x and dbinom(x = 10000, size = 1, prob = 0.4)
## D = 0.3941, p-value = 0.9977
## alternative hypothesis: two-sided

Punto 4

**Solucíon:

para dar solución, pensemos en el metodo de la inversión.

rPoisson <- function(lambda = 1){
  U <- runif(1)
  i <- 0
  p <- exp(-lambda)
  f <- p
  while(U >= f){
    p <- lambda * p / (i + 1)
    f <- f + p
    i <- i + 1
  }
  i
}
sim5000<- replicate(5000,rPoisson())
sim10000<-replicate(1000,rPoisson())

plot(factor(sim5000))

plot(factor(sim10000))

Veamos la grafica para los numeros aleatorios “teoricos” generados con R, y veamos que gráficamente son igules.

x<-rpois(5000,1)
plot(factor(x))

x<-rpois(10000,1)
plot(factor(x))

Prueba de bondad de ajuste: para n=5000 y n=10000

ks.test(sim5000,dpois(5000,1))
## Warning in ks.test(sim5000, dpois(5000, 1)): cannot compute exact p-value with
## ties
## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  sim5000 and dpois(5000, 1)
## D = 0.6292, p-value = 0.8235
## alternative hypothesis: two-sided
ks.test(sim10000,dpois(10000,1)) 
## Warning in ks.test(sim10000, dpois(10000, 1)): cannot compute exact p-value with
## ties
## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  sim10000 and dpois(10000, 1)
## D = 0.637, p-value = 0.8123
## alternative hypothesis: two-sided

Luego los datos generados si se ajustan a una distribución Poisson, por que el p>0.05