TALLER 2 DE SIMULACIÓN

PUNTO 1

Dado el siguiente algoritmo:

Genere un numero aleatorio \(u\)~\(unif(0,1)\)
Si \(u < p_{0}\) hacer \(x = x0\)
Si no, si \(u < p_{0} + p_{1}\) hacer \(x = x_{1}\)
Si no, si \(u < p_{0} + p_{1} + p_{2}\) hacer \(x = x_{2} ...\)
Si no, si \(u < p_{0} + p_{1} + p_{2} + + p_{j}\) hacer \(x = x_{j}\)
Fin del siclo
Volver al paso \(1\) \(n-veces.\)

Simular para dos tamaños de muestra \(n=5000\) y \(n=10000\) variables aleatorias de tal modo que:

\(p_{1}=[x=1]=0.20\)

\(p_{1}=[x=1]=0.15\)

\(p_{1}=[x=1]=0.25\)

\(p_{1}=[x=1]=0.40\)

Para cada tamaaño de muestra realizar un diagrama de barras.

Solucion:

Para muesta de tamaño n=5000

n=5000
p1=0.20
p2=0.15
p3=0.25
p4=0.40
X=1

for (i in 1:n) {
  
  U<-runif(1)
  if(U<p1){X[i]=1
  }else{
    if(U>=p1 & U<p1+p2){
      X[i]=2
    }else{
      if(U>=p1+p2 & U<p1+p2+p3){
        X[i]=3
      }else{
        if(U>=p1+p2+p3 & U<p1+p2+p3+p4){
          X[i]=4
        }
      }
    }
  }
  X
}
#X
table(X)

## X
##    1    2    3    4 
##  993  767 1238 2002

plot(factor(X),main='Diagrama de barras',col=38)

Para muestra de tamaño n=10000

n=10000
p1=0.20
p2=0.15
p3=0.25
p4=0.40
X=1

for (i in 1:n) {
  
  U<-runif(1)
  if(U<p1){X[i]=1
  }else{
    if(U>=p1 & U<p1+p2){
      X[i]=2
    }else{
      if(U>=p1+p2 & U<p1+p2+p3){
        X[i]=3
      }else{
        if(U>=p1+p2+p3 & U<p1+p2+p3+p4){
          X[i]=4
        }
      }
    }
  }
  X
}
#X
table(X)

## X
##    1    2    3    4 
## 2077 1451 2540 3932

plot(factor(X),main='Diagrama de barras',col=34)

PUNTO 2

Dado el siguiente algoritmo:

Hacer \(j =\) valor maximo de la variable aleatoria.
Generar un numero aleatorio u~unif(0,1)
Hacer \(x = ent(j*u) + 1\)
Volver al paso 2 n veces

Simular para dos tamaños de muestra \(n=5000\) y \(n=10000\) variables aleatorias que se comporten como una distribucion uniforme discreta, convertir el programa en una funcion. Para cada tamaño de muestra realizar un diagrama de barras y una prueba de bondad de ajuste.

Solucion:

Para n=5000:

set.seed(5)
z<-1
Unifdis<-function(n,j){
n<-5000
for (i in 1:n) {
x<-c(1,2,3,4,5,6)
j<-max(x)
u<-runif(1)
z[i]<-as.integer(j*u)+1
}
return(z)
}
x<-Unifdis(100,6)
table(x)

## x
##   1   2   3   4   5   6 
## 794 845 851 829 822 859

plot(factor(x),main='Diagrama de barras',col=34)

Prueba de bondan de ajuste:

chisq.test(table(x))

## 
##  Chi-squared test for given probabilities
## 
## data:  table(x)
## X-squared = 3.3616, df = 5, p-value = 0.6444

Conclusion:

Como el valor de p = 0.6444 es mayor que el de \(\alpha= 0.05\) no se rechaza \(H_{o}\). Los datos de la muestra se ajustan a una distribución uniforme discreta.
Como el \(3.3616\) no es mayor a \(11.07\), no se rechaza H0. y se concluye con un \(p = 0.05\) que el ajuste de los datos a una distribución uniforme discreta es bueno.

Para n=10000:

set.seed(8)
z<-1
Unifdis<-function(n,j){
n<-10000
for (i in 1:n) {
x<-c(1,2,3,4,5,6)
j<-max(x)
u<-runif(1)
z[i]<-as.integer(j*u)+1
}
return(z)
}
x<-Unifdis(100,6)
table(x) #tabla de frecuencia

## x
##    1    2    3    4    5    6 
## 1679 1677 1627 1652 1659 1706

plot(factor(x),main='Diagrama de barras',col=69)

Prueba de bondan de ajuste:

chisq.test(table(x))

## 
##  Chi-squared test for given probabilities
## 
## data:  table(x)
## X-squared = 2.192, df = 5, p-value = 0.822

Conclusion:

Como el valor de p = 0.822 es mayor que el de p = 0.05 no se rechaza Ho. Los datos de la muestra se ajustan a una distribución uniforme discreta.
Como el \(2.192\) no es mayor a \(11.07\), no se rechaza \(H_{0}\). y se concluye con un \(\alpha = 0.05\) que el ajuste de los datos a una distribución uniforme discreta es bueno.

PUNTO 3

Dado el siguiente algoritmo que genera variables aleatorias por el metodo de Aceptacion y Rechazo:

Simular un valor y que tenga función de probabilidad \(q\).
Genere un número aleatorio u.
Si \(u <= (p_{y}/c*q_{y})\) asignar \(x = y\)
Si no, volver al paso 1

Simular y graficar para dos tamaaños de muestra \(n=5000\) y \(n=10000\) variables aleatorias de tal modo que:

\(p_{1} = [x = 1] = 0,11\)

\(p_{2} = [x = 2] = 0,12\)

\(p_{3} = [x = 3] = 0,09\)

\(p_{4} = [x = 4] = 0,08\)

\(p_{5} = [x = 5] = 0,12\)

\(p_{6} = [x = 6] = 0,10\)

\(p_{7} = [x = 7] = 0,09\)

\(p_{8} = [x = 8] = 0,09\)

\(p_{9} = [x = 9] = 0,10\)

\(p_{10}= [x = 10]= 0,10\)

Solucion:

Para n=5000

set.seed(1)
n<-5000
x <- c(1,2,3,4,5,6,7,8,9,10)
lx<-length(x)
p <- c(0.11, 0.12,0.09,0.08,0.12,0.10,0.09,0.09,0.10,0.10)
q <- rep(1/lx, lx)
Max<-max(p)/max(q)
X <- rep(0, n)

for (i in 1:n) {
  for (j in 1:100) {
    Y <- sample(1:lx, 1)
    U <- runif(1)
    if (U < p[Y]/((Max) * q[Y])) {
      X[i] <- x[Y]
    }
  }
}
#X
table(X)

## X
##   1   2   3   4   5   6   7   8   9  10 
## 557 592 482 389 624 485 454 432 496 489

plot(factor(X),main='Diagrama de barras',col=34)

Para n=10000

set.seed(1)
n<-10000
x <- c(1,2,3,4,5,6,7,8,9,10)
lx<-length(x)
p <- c(0.11, 0.12,0.09,0.08,0.12,0.10,0.09,0.09,0.10,0.10)
q <- rep(1/lx, lx)
X <- rep(0, n)

for (i in 1:n) {
  for (j in 1:100) {
    Y <- sample(1:lx, 1)
    U <- runif(1)
    if (U < p[Y]/(2 * q[Y])) {
      X[i] <- x[Y]
    }
  }
}
#X
table(X)

## X
##    1    2    3    4    5    6    7    8    9   10 
## 1084 1160  915  787 1227 1014  909  910  962 1032

plot(factor(X),main='Diagrama de barras',col=15)

PUNTO 4

Usar simulacion para aproximar las siguientes integrales para tres tamaños de muestras. Compare su estimacion con el valor exacto.

\(a.\int_{0}^{1}e^{-{x/2}}dx\)

*Hallamos el valor de la integral aplicando los comandos de R.

Int1<-function(x){(exp(-(x/2)))}

integrate(Int1,0,1) #la integral

## 0.7869387 with absolute error < 8.7e-15

Hallamos el valor de la integral por Simulación Montecarlo:

set.seed(5)
a=0
b=1
nsim=10000000
U=runif(nsim,0,1)
I=(b-a)*(exp(-((U*(b-a)+a)/2)))

mean(I) ####corresponde al valor estimado de la integral

## [1] 0.7869912

hist(I,col=75)

\(b.\int_{2}^{4}e^{-x/}dx\)

*Hallamos el valor de la integral aplicando los comandos de R.

Int2<-function(x) {exp(-x)}

integrate(Int2,2,4) ##corresponde al valor teorico de la integral

## 0.1170196 with absolute error < 1.3e-15

curve(Int2,xlim=c(-4,4),lwd=3)

Hallamos el valor de la integarl por Simulacion Montecarlo:

set.seed(5)

a<-2
b=4
nsim=10000000
U=runif(nsim,0,1)
I2=(b-a)*(exp(-(U*(b-a)+a)))

mean(I2) ##corresponde al valor estimado de la integral

## [1] 0.1170526

hist(I2,col=75)

Hallamos el valor de la integarl por Simulacion Montecarlo Por una segunda forma:

set.seed(5)
a=2
b=4
nsim=100000 
x=matrix(0,1,nsim)
 for (i in 1:nsim){
u=runif(1)
x[1,i]=((b-a)*(exp(-(u*(b-a)+a))))
 }

mean(x) ##corresponde al valor estimado de la integral

## [1] 0.11758

hist(x, col= 76)

\(c.\int_{0}^{\infty}e^{-2x}dx\)

Hallamos el valor de la integral aplicando los comandos de R.

Int3<-function(x){exp(-2*x)}

integrate(Int3,0,Inf) #Corresponde al valor teorico de la integral

## 0.5 with absolute error < 7.7e-05

Hallamos el valor de la integarl por Simulacion Montecarlo:

a=-7
b=-5
nsim=10000000
U=runif(nsim,0,1)
I3=(exp(-2*((1/U)-1))/U^2)

mean(I3)##corresponde al valor estimado de la integral

## [1] 0.4998591

hist(I3,col=75)

\(d.\int_{0}^{1}\int_{0}^{1}(2x+6x^2y)dydx\)

Hallamos la integral por similacion Montecarlo:

set.seed(7)
nsim=100000
u1=runif(nsim,0,1)
u2=runif(nsim,0,1)
Int4=2*u1+6*(u1^2)*u2
mean(Int4)  ##corresponde al valor estimado de la integral

## [1] 2.000845

hist(Int4,col=67, main='')

PUNTO 5

Simular el teorema del límite central para dos distribuciones una discreta y una continua. 1. Generar una muestra aleatoria con distribución Exponencial \((n,p)\) de tamaño \(n1\).

Calcular \(\bar{X}_{n1}\)
\(Z= (\bar{X}_{n1}-\mu )/(\sigma /\sqrt{n1})\)
Guardar \(Z\) en \(Y\)
Volver al paso 1 \(n2-veces\)

Simular para diferentes valores de \(n1\) y \(n2\), realizar gráficas y pruebas de bondad de ajuste.

Solucion:

Distribucion exponencial

set.seed(1)
TLC_Exp=function(N,a,n){
  mu=1/a
  var=1/a^2
  y=c()
  for(i in 1:N){
    s=sqrt(var)
    Exp=rexp(N,a)
    m=mean(Exp)
    y[i]=((m-mu)/(s/sqrt(N)))
  }
  y
}

X1<-TLC_Exp(10000,1,30)
X2<-TLC_Exp(2000,1,20)
X3<-TLC_Exp(6000,1,40)

par(mfrow = c(2,2))
hist(X1,col=76) ##Para observar que la media se distribuye como una normal 
hist(X2,col=7)
hist(X3,col=27)

ks.test(X,dexp(1000,1)) #analogamente se hace la prueba para X2 y X3

## Warning in ks.test(X, dexp(1000, 1)): p-value will be approximate in the
## presence of ties

## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  X and dexp(1000, 1)
## D = 1, p-value = 0.2701
## alternative hypothesis: two-sided

Conclusión:

Como el valor de \(p = 0.96\) es mayor que el de \(\alpha = 0.05\) no se rechaza Ho. Los datos de las muestras de las medias se ajustan a una distribución normal.

Punto 6

Simular la desigualdad de Chebyshev siguiendo el siguiente algoritmo:

Generar una muestra aleatoria de tamaño \(n\) de una normal de media \(\mu\) y desviación estándar \(\sigma\).
Calcular \(\bar{X}_{n}\), hacer \(suma=0\)
Hallar \(\bar{X}_{n}-\mu<=\delta\) \(suma=suma+1\)
Si no, vuelvo al paso \(2\)
Repito \(N-veces\).
Calculo probabilidad, \(suma/N\).

y<-0
nsim<-1000

for (i in 1:nsim) {
  k<-2
  n<-1 #tamaño de muestra
  m<-50
  sd<-5
  Delta<- k*sd
x<-rnorm(n,m,sd)
xn<-mean(x)

if(abs(xn-m)<=Delta)
{y[i]=1

  }else{if(abs(xn-m)>Delta){y[i]=0}}

}
#y
table(y)

## y
##   0   1 
##  47 953

sum(y)

## [1] 953

sum(y)/nsim #Esta es laprobabilidad pedida

## [1] 0.953

TALLER 2 DE SIMULACIÓN

Jose David Ortega B- (070250112017) y Antenor Amezquita Orrego-(0700200192015)

19/8/2022

PUNTO 1

Solucion:

PUNTO 2

Solucion:

Para n=5000:

Para n=10000:

PUNTO 3

Solucion:

Para n=5000

Para n=10000

PUNTO 4

\(a.\int_{0}^{1}e^{-{x/2}}dx\)

\(b.\int_{2}^{4}e^{-x/}dx\)

\(c.\int_{0}^{\infty}e^{-2x}dx\)

\(d.\int_{0}^{1}\int_{0}^{1}(2x+6x^2y)dydx\)

PUNTO 5

Solucion:

Distribucion exponencial

Punto 6