Ejemplo 1

n_caras = 6; n_lanzamientos = 1200
resultado = c()

for (i in 1:n_lanzamientos) {
  resultado[i] = sample(1:n_caras, 1)
}
#Imprimiento los resultados 
resultado
##    [1] 1 1 6 5 2 5 1 5 6 5 6 4 2 5 4 4 1 6 6 2 2 3 1 3 1 4 1 5 1 2 4 5 4 1 2 6 3
##   [38] 1 1 3 6 6 3 4 4 5 2 4 6 1 1 4 6 3 1 4 2 1 3 3 6 3 2 5 3 6 6 5 2 6 6 6 3 1
##   [75] 2 2 4 6 1 1 6 4 3 4 4 5 1 6 3 2 6 5 3 5 3 6 3 1 5 2 6 4 2 6 2 3 5 2 5 2 2
##  [112] 4 6 3 4 5 6 3 5 4 6 1 6 6 1 1 2 5 1 1 1 2 6 6 1 5 4 4 6 6 5 1 1 4 3 3 3 5
##  [149] 6 1 3 6 4 4 6 3 5 5 5 3 6 6 2 2 5 2 3 3 1 3 4 3 2 3 1 1 1 6 2 4 1 1 5 5 2
##  [186] 3 1 5 1 5 5 2 4 1 3 2 1 1 3 6 2 3 5 6 3 6 5 1 1 5 4 6 6 2 6 4 4 5 1 2 4 4
##  [223] 1 4 2 6 1 5 6 4 1 3 5 1 3 2 2 4 1 3 5 2 2 3 6 5 1 2 4 1 3 3 4 3 3 3 5 3 4
##  [260] 3 3 6 2 2 2 5 5 6 4 1 5 3 6 1 4 4 6 4 3 4 6 4 3 3 2 2 2 6 2 4 6 1 6 5 2 4
##  [297] 5 2 4 6 5 6 2 5 2 5 1 4 5 3 1 1 2 3 4 5 6 5 5 3 2 4 1 1 1 6 6 1 1 1 3 2 1
##  [334] 5 5 2 1 4 1 5 4 1 2 1 3 2 2 4 3 2 3 3 2 1 3 1 4 4 1 6 3 4 3 3 5 3 6 5 6 6
##  [371] 5 4 3 3 4 2 6 5 6 6 4 3 2 4 6 3 1 3 3 4 1 2 2 2 1 2 2 1 3 6 5 5 6 2 4 2 4
##  [408] 6 1 4 2 6 5 6 5 5 6 3 1 6 5 1 4 5 5 6 4 5 2 6 6 3 6 4 2 5 6 6 6 1 3 6 4 3
##  [445] 2 3 3 3 5 4 5 5 6 5 5 1 6 2 5 5 2 3 4 5 3 4 6 4 3 6 5 1 4 6 3 5 3 5 4 3 2
##  [482] 1 2 6 4 1 3 5 2 6 6 3 1 1 5 6 5 5 2 2 3 2 6 2 3 5 2 2 4 5 1 3 1 3 5 3 3 1
##  [519] 5 6 3 2 5 2 4 4 2 5 3 5 4 5 3 3 5 6 3 6 4 2 5 6 4 4 1 1 4 6 5 3 3 4 2 6 4
##  [556] 6 6 3 1 5 1 4 6 2 6 1 4 4 5 4 3 3 6 1 3 3 2 2 2 5 5 2 5 6 4 1 5 6 5 2 6 5
##  [593] 6 1 6 3 3 3 4 3 5 1 2 3 2 1 4 1 2 6 2 5 2 2 3 5 3 6 3 3 4 5 5 6 1 4 5 6 2
##  [630] 3 4 2 6 1 4 2 5 5 1 5 2 2 3 5 5 4 6 4 1 1 1 6 1 5 3 5 6 5 5 5 3 3 3 3 3 1
##  [667] 5 6 3 1 2 5 4 6 1 6 4 3 4 3 3 1 2 1 4 3 2 5 4 6 5 1 4 4 1 6 5 1 2 3 3 3 4
##  [704] 2 3 2 4 5 1 1 1 4 1 4 6 2 6 5 3 6 3 2 2 6 2 2 4 5 5 2 1 3 3 6 5 2 4 3 1 4
##  [741] 5 2 6 6 5 3 5 2 5 4 5 5 3 5 5 3 3 5 2 4 1 4 2 1 1 2 5 3 3 3 5 1 3 4 5 4 4
##  [778] 5 6 5 5 4 5 2 2 6 5 3 5 4 6 2 3 2 3 3 2 1 2 2 2 6 3 4 4 3 1 6 5 3 2 2 2 6
##  [815] 2 5 6 5 2 3 6 2 2 3 3 5 3 2 5 4 6 5 5 2 6 6 4 4 5 6 2 1 4 2 4 6 4 1 6 6 4
##  [852] 1 6 5 6 6 3 1 3 6 1 3 2 2 2 5 2 3 3 3 6 6 3 2 4 6 6 3 4 4 3 3 2 2 5 3 5 3
##  [889] 6 2 6 6 2 5 2 1 1 4 3 4 2 5 3 4 3 2 3 2 3 5 5 4 6 3 3 2 6 2 6 2 5 6 4 3 5
##  [926] 4 4 6 1 4 1 6 1 3 1 2 3 6 3 6 6 2 4 5 4 6 1 2 1 1 3 2 1 1 1 3 3 2 2 6 2 2
##  [963] 3 2 2 4 2 4 6 6 3 2 1 6 2 2 6 6 6 6 6 2 4 2 1 6 2 3 3 5 4 1 4 5 2 6 4 6 3
## [1000] 1 4 1 5 4 3 1 6 5 5 5 4 6 2 5 3 6 3 4 6 5 4 1 3 2 6 3 6 3 3 6 3 6 2 3 2 1
## [1037] 5 4 6 5 4 1 3 1 6 1 4 4 5 2 4 1 6 2 4 3 6 4 2 2 3 5 1 4 6 5 3 3 6 1 4 2 1
## [1074] 2 6 2 5 4 4 1 6 2 2 1 5 3 3 4 2 5 5 3 3 2 5 2 6 6 1 3 6 5 5 1 5 5 2 2 2 3
## [1111] 3 3 1 4 5 4 5 5 5 2 3 2 5 6 3 3 5 6 5 6 2 5 2 6 6 4 1 6 5 3 2 1 2 4 3 5 5
## [1148] 5 5 2 3 5 3 3 4 1 3 4 5 6 1 6 3 1 5 4 4 6 6 6 4 4 5 2 1 3 1 1 6 2 1 4 1 5
## [1185] 5 1 4 2 2 3 4 2 6 6 4 1 6 1 4 1
#Tabulando los resultados
table(resultado)
## resultado
##   1   2   3   4   5   6 
## 177 206 219 180 207 211
#Proporcion de datos
prop.table(table(resultado))
## resultado
##         1         2         3         4         5         6 
## 0.1475000 0.1716667 0.1825000 0.1500000 0.1725000 0.1758333

Ejemplo 2

n_caras = 6; n_lanzamientos = 12000
resultado_1 = c()
resultado_2 = c()
suma = c()
prod = c()

for (i in 1:n_lanzamientos) {
  resultado_1[i] = sample(1:n_caras, 1)
  resultado_2[i] = sample(1:n_caras, 1)
  suma[i] = resultado_1[i] + resultado_2[i]
  prod[i] = resultado_1[i] * resultado_2[i]
}

#Tabulando los resultados
tabla1= table(resultado_1)
tabla2= table(resultado_2)
#Proporcion de datos
prop.table(tabla1)
## resultado_1
##         1         2         3         4         5         6 
## 0.1700833 0.1682500 0.1605000 0.1699167 0.1646667 0.1665833
prop.table(tabla2)
## resultado_2
##         1         2         3         4         5         6 
## 0.1637500 0.1644167 0.1677500 0.1684167 0.1706667 0.1650000
#Suma mas probable
barplot(table(suma))

#Producto mas probable
barplot(table(prod))

#Tabla cruzada
tabla.cruzada = table(resultado_1,resultado_2)
#Porcentaje de celda
prop.table(tabla.cruzada)
##            resultado_2
## resultado_1          1          2          3          4          5          6
##           1 0.02483333 0.02866667 0.03125000 0.03050000 0.02650000 0.02833333
##           2 0.02675000 0.02766667 0.02750000 0.02900000 0.02625000 0.03108333
##           3 0.02808333 0.02516667 0.02566667 0.02558333 0.02991667 0.02608333
##           4 0.02991667 0.02816667 0.02825000 0.02850000 0.02658333 0.02850000
##           5 0.02683333 0.02841667 0.02650000 0.02750000 0.03125000 0.02416667
##           6 0.02733333 0.02633333 0.02858333 0.02733333 0.03016667 0.02683333
#Porcentaje de columna
prop.table(tabla.cruzada,2)
##            resultado_2
## resultado_1         1         2         3         4         5         6
##           1 0.1516539 0.1743538 0.1862891 0.1810985 0.1552734 0.1717172
##           2 0.1633588 0.1682717 0.1639344 0.1721920 0.1538086 0.1883838
##           3 0.1715013 0.1530664 0.1530055 0.1519050 0.1752930 0.1580808
##           4 0.1826972 0.1713127 0.1684054 0.1692232 0.1557617 0.1727273
##           5 0.1638677 0.1728332 0.1579732 0.1632855 0.1831055 0.1464646
##           6 0.1669211 0.1601622 0.1703924 0.1622959 0.1767578 0.1626263

Ejemplo 3

n = 10000
x = runif(n, 0, 1)
y = runif(n, 0, 1)

d = (x^2+y^2<=1)
plot(x,y, pch= 16, asp = TRUE, col=d, cex=0.1)
abline(v=c(0,1),h=c(0,1), col='green')

#Se tabula "d" para saber cuantos puntos caen dentro del circulo
table(d)
## d
## FALSE  TRUE 
##  2235  7765
#Porcentaje de los puntos que caen dentro
100*table(d)/n
## d
## FALSE  TRUE 
## 22.35 77.65

Si el area completa del cuadrado es 1 u^2 entonces 79% del area son los puntos <= 1, es decir 0.79 unidades de area es el area de ese cuarto de circulo

#Calculando el area real 
a_real = (pi*1^2)/4;a_real
## [1] 0.7853982
#Calculando area estimada a traves de simulacion
a_esti = sum(d)/n;a_esti
## [1] 0.7765

Ejemplo 4

n = 10000
cociente = c()
for(i in 10:n){
  x = runif(i, 0, 1)
  y = runif(i, 0, 1)
  
  d = (x^2+y^2<=1)
  area_esti = sum(d)/i
  area_real = pi/4
  
  cociente[i-9] = area_esti/area_real
}
length(cociente)
## [1] 9991
#Cociente en la posicion 2000
cociente[2000]
## [1] 1.013395
#Cociente en la posicion 5000: no se requiere un numero de simulaciones tan grande para que el resultado sea cercano a 1
cociente[5000]
## [1] 0.9964262
#Cociente en desde la posicion 20 hasta la 30
cociente[20:30]
##  [1] 0.8780962 0.9761503 1.1500228 0.9549297 1.0803245 1.0859984 1.0549699
##  [8] 0.9902974 0.9291207 0.9381765 1.1100037
#Graficando el cociente 
plot(cociente, pch=16, cex=0.1)
abline(h=1, col='green')

Asignacion 1 Estimar el area de una elipse cuyo (semieje mayor) a = 2 y (semieje menor) b = 1

n = 100000
x= runif(n, -2, 2)
y= runif(n, -1, 1)
e = ((x^2/4+y^2)<=1)
plot(x,y, pch= 16, asp = TRUE, col=e, cex=0.1)
abline(v=c(-2,2),h=c(-1,1), col='green')

#Se tabula "d" para saber cuantos puntos caen dentro del circulo
table(e)
## e
## FALSE  TRUE 
## 21486 78514
#Porcentaje de los puntos que caen dentro
100*table(e)/n
## e
##  FALSE   TRUE 
## 21.486 78.514
#Calculando el area real 
a_real = (pi*2*1);a_real
## [1] 6.283185
#Calculando area estimada a traves de simulacion
a_esti = sum(e)/n;a_esti
## [1] 0.78514

Ejemplo 5: area de una integral por simulacion

curve(sqrt(x), from = 0, to = 10)
abline(v = c(0,10), h=0, col=gray(0.7))
segments(6,0,6,sqrt(6), col='red', lty=2) 
segments(6.5,0,6.5,sqrt(6.5), col='red', lty=2)
text(6.3,2.7, 'dx')

\[A=\int_{0}^{10}f(x)dx\] \[A=\int_{0}^{10}\sqrt{x}dx\] \[A=\left [\frac{2}{3} x^{\tfrac{3}{2}} \right ]_{0}^{10}\] \[A=\frac{2}{3} 10^{\tfrac{3}{2}}\] \[A=21.082\]

#Estimacion para hallar el area
n=100000
x = runif(n, 0, 10)
y = runif(n, 0, sqrt(10))

punt = (y<=sqrt(x))
plot(x,y, pch=16, cex=0.1, col=punt+2)
curve(sqrt(x), add=T, col='blue')

#Area total de todo el cuadrado 
area_cuadrado = 10*sqrt(10); area_cuadrado
## [1] 31.62278
#El 66% de los puntos cayeron debajo de la curva
puntos_curva =100*(sum(punt)/n);puntos_curva
## [1] 66.579
#Area estimada 
area_estimada =(sum(punt)/n)*area_cuadrado;area_estimada
## [1] 21.05413

Asignacion 2 Hallar el area encerrada entre las curvas siguiente y compararla con la real calcula por integrales

curve(expr = sqrt(x), from = 0, to = 1.5, ann=F)
curve(expr = x**2, from = 0, to = 1.5, add=T)
text(0.25,0.25,'A')
text(0.25,0.8,expression(sqrt(x)))
text(0.7,0.25,expression(x^2))

#Calcular integral
integrate(function(x)sqrt(x)- x^2, lower = 0, upper = 1)
## 0.3333334 with absolute error < 7.8e-05
#Estimacion para hallar el area entre curvas
n=100000
x = runif(n, 0, 1)
y = runif(n, 0,1)


puntos= (y<=sqrt(x) & y>=(x**2))
plot(x,y, pch=16, cex=0.1,col=puntos)
curve(sqrt(x), add=T, col='red')
curve(x**2, add=T, col='red')

#Area total de todo el cuadrado 
area_cuadrado = 1*1; area_cuadrado
## [1] 1
#El 33% de los puntos cayeron debajo de la curva
puntos_curva =100*(sum(puntos)/n);puntos_curva
## [1] 33.213
#Area estimada 
area_estimada =(sum(puntos)/n)*area_cuadrado;area_estimada
## [1] 0.33213