Estimación del valor de π. LA siguiente figura sugiere como estimar el valor de πcon una simulación. En la figura, un circuíto con un áreaigual a π/4, está inscrito en un cuadrado cuya área es igual a 1. Se elige de forma aleatoria 100 puntos dentro del cuadrado . La probabilidad de que un punto esté dentro del círculo es igual a la pracción del área del cuadrado que abarca a este, la cual es π/4. Por tanto, se puede estimar el valor de π/4 al contar el número de puntos dentro del círculo, que es 79 para obtener la estimación de π/4≈0.76. De este último resultado se concluye que π≈4(0.79)=3.14. Este ejercicio presenta un experimento de simulación que fue diseñado para estimar el valor de π al generar 1000 puntos en el cuadrado.
Genere 1000 coordenadas x: X1, . . . , X1000. Utilice la distribución uniforme con valor mínimo de 0 y valor máximo de 1. La distribución uniforme genera variables aleatorias que tienen la misma probabilidad de venir de cualquier parte del intervalo (0,1).
x = runif(1000,0,1)
Genere 1000 coordenadas y : Y1,…,Y1000, utilizando nuevamente la distribución uniforme con valor mínimo de 0 y valor máximo de 1.
y = runif(1000,0,1)
Cada punto (Xi,Yi) se encuentra dentro del círculo si su distancia desde el centro (0.5,0.5) es menor a 0.5. Para cada par (Xi,Yi) determine si la distancia desde el centro es menor a 0.5. Esto último se puede realizar al calcular el valor (Xi−0.5)^2 + (Yi−0.5)^2, que es el cuadrado de la distancia, y al determinar si es menor que 0.25.
#creamos un vector que almacena las distancias
dis_p_c <- numeric()
#calculamos las distancias
for (i in 1:length(x)){
dis_p_c[i] = (x[i]-0.5)^2 + (y[i]-0.5)^2
}
#guardamos los puntos tanto dentro como fuera del circulo
#1 -> esta dentro
#0 -> esta fuera
puntos = ifelse(dis_p_c < 0.25, 1, 0)
¿Cuántos de los puntos están dentro del círculo? ¿Cuál es su estimación de π? (Nota: Con sólo 1000 puntos, es probable que su estimación sea inferior por 0.05 o más. Una simulación con 10000 y 100000 puntos tiene mayores probabilidades de dar como resultado una estimación muy cercana al valor verdadero
#graficamos los puntos
t <- seq(0, 2*pi, length.out = 100)
plot(0, 0, asp = 1, type = "n",
xlim = c(0, 1), ylim = c(0, 1),
ann = F)
abline(v = seq(0, 1, 0.1), lty = 2, col = "gray50")
abline(h = seq(0, 1, 0.1), lty = 2, col = "gray50")
radio <- 1/2
a <- 0.5 ## origen circunferencia eje x
b <- 0.5 ## origen circunferencia eje y
xx <- a + cos(t)*radio
yy <- b + sin(t)*radio
points(xx, yy, type = "l", col = "blue")
points(x,y,pch = 20)
#numero de puntos dentro del circulo: sumamos cuantos valores de puntos son iguales a 1.
sum(puntos)
## [1] 788
#estimamos pi
4 * sum(puntos)/length(x)
## [1] 3.152
4 * sum(puntos)/length(x) - pi
## [1] 0.01040735
La edad de una antigua pieza de materia orgánica se puede estimar a partir de la tasa a la que emite partículas beta como resultado del decaimiento del carbono-14. Por ejemplo, si X es el número de partículas emitidas durante diez minutos por un fragmento óseo con 10000 años de antigüedad que contiene 1 g de carbono, entonces X tiene una distribución de Poisson con media λ=45.62 . Un arqueólogo descubrió un pequeño fragmento óseo que contiene 1 g de carbono. Si t^es la edad desconocida del hueso, en años, el arqueólogo contará el número X de partículas emitidas en diez minutos y calculará una edad estimada t^con la fórmula:
El arqueólogo no lo sabe,
pero el hueso tiene exactamente 10000 años de antiguedad, por lo que X
tiene una distribución de Poisson con λ=45.62.
a.Genere una muestra simulada de 10000 valores de X y sus correspondientes valores de t^.
X = rpois(10000, 45.62)
t = (log(15.3) - log(X/10)) / 0.0001210
d = data.frame(
X = X,
t = t
)
View(d)
b.Estime la media de t^.
media = mean(t)
media
## [1] 10107.99
c.Estime la desviación estándar de t^.
sqrt(media)
## [1] 100.5385
d.Estime la probabilidad de que t^esté a 1000 años con una edad real de 10000 años.
lower_bound <- 10000 - 1000
upper_bound <- 10000 + 1000
prob <- mean(t >= lower_bound & t <= upper_bound)