GeneraciĆ³n de Variables Aleatorias no Uniforme
Odacir Padilla S.
2024-09-18
IntroducciĆ³n
En el proceso de realizar una simulaciĆ³n de eventos necesitaremos conocer la distribuciĆ³n de probabilidad que modela dicho evento , por ejemplo si nuestro propĆ³sito al analizar un muestreo de calidad consiste en modelar \(X\) nĆŗmero de artĆculos defectuosos en una muestra de \(n\) artĆculos, entonces \(X\sim Binomial(n,p)\) , donde \(p\) es la probabilidad de que el artĆculo este defectuoso.
En cambio si , estamos analizando \(Y:\) nĆŗmero de usuarios que llamarĆ”n a un telĆ©fono de atenciĆ³n a clientes en un determinado tiempo,entonces \(X\sim Poisson(\lambda)\), y si estamos interasados analizar la variables \(T:\) _Tiempo entre llegada de clientes__, entonces la variable \(T\sim exp(\beta)\),
a estas distribuciĆ³n de variables aleatorias se asocia una funciĆ³n \(F\) llamada funciĆ³n de distribuciĆ³n acumulada \[F(x)=P(X\leq x)\] Para simular variables aleatorias \(X\) de una determinada distribuciĆ³n \(f(x)\) , primero necesitamos un generador de nĆŗmeros pseudoaleatorios que produzca valores uniformemente distribuidos. Luego, mediante un mĆ©todo especĆfico, transformamos esos nĆŗmeros uniformes en valores que sigan la distribuciĆ³n deseada.
- MĆ©todo de la Tranformada Inversa
Es una tĆ©cnica para generar variables aleatorias de cualquier distribuciĆ³n de probabilidad a partir de nĆŗmeros aleatorios uniformemente distribuidos. Es especialmente Ćŗtil cuando la funciĆ³n de distribuciĆ³n acumulada \(F\) de la variable deseada es invertible, es decir existe \(F^{-1}\)
Pasos:
Se genera un nĆŗmero aleatorio \(U\) a partir de una distribuciĆ³n uniforme en el intervalo [0,1].
Se utiliza la funciĆ³n de distribuciĆ³n acumulada \(š¹(š„)\) de la variable aleatoria deseada.
Se obtiene \(F^{-1}(u)\) que permite transformar un nĆŗmero uniforme \(U\) en un valor que sigue la distribuciĆ³n deseada.
El nĆŗmero uniforme \(U\) se transforma en una variable aleatoria \(X\) utilizando la inversa \(F^{-1}\), es decir: \[X=F^{-1}(U)\] Este \(X\) sigue la distribuciĆ³n de probabilidad deseada.
Ejemplo.
Para una variablea aleatoria \(X\sim exp(\beta)\) su funciĆ³n de distribuciĆ³n acumulada \(F\) es :
\[F(x)= 1-exp(-\frac{1}{\beta}x)\]
Para aplicar el mƩtodo de la transformada inversa:
- Generamos un nĆŗmero \(Uā¼U(0,1)\).
- Hacemos \(U=1-exp(-\frac{1}{\beta}x)\)
- despejamos \(x\)
\[x=-\beta \ln(1-u)\] \(x\) sigue una distribuciĆ³n exponenial con parĆ”metro \(\beta\).
Los datos histĆ³ricos del tiempo de servicio en la caja de un banco se com portan de forma exponencial con media de 3 minutos/cliente. con los nĆŗmeros aleatorios generados por el mĆ©todo Congriuencial mixto.
a <- 317
c <- 15
m <- 571
X_n <- 41 # semilla
random.number<-numeric(50) # vector numƩrico de longitud 50
for (i in 1:50)
{X_n<-(a*X_n+c)%%m
random.number[i]<-X_n/m # nĆŗmeros en el intervalo [0,1]
}
random.number## [1] 0.78809107 0.85113835 0.83712785 0.39579685 0.49387040 0.58318739
## [7] 0.89667250 0.27145359 0.07705779 0.45359019 0.81436077 0.17863398
## [13] 0.65323993 0.10332750 0.78108581 0.63047285 0.88616462 0.94045534
## [19] 0.15061296 0.77057793 0.29947461 0.95971979 0.25744308 0.63572680
## [25] 0.55166375 0.90367776 0.49211909 0.02802102 0.90893170 0.15761821
## [31] 0.99124343 0.25043783 0.41506130 0.60070053 0.44833625 0.14886165
## [37] 0.21541156 0.31173380 0.84588441 0.17162872 0.43257443 0.15236427
## [43] 0.32574431 0.28721541 0.07355517 0.34325744 0.83887916 0.95096322
## [49] 0.48161121 0.69702277Entonces
# generando una variable exponencial con beta = 3
x <- -3*log(1-random.number)
round(x,1) # redondeando a una cifra decimal## [1] 4.7 5.7 5.4 1.5 2.0 2.6 6.8 1.0 0.2 1.8 5.1 0.6 3.2 0.3 4.6
## [16] 3.0 6.5 8.5 0.5 4.4 1.1 9.6 0.9 3.0 2.4 7.0 2.0 0.1 7.2 0.5
## [31] 14.2 0.9 1.6 2.8 1.8 0.5 0.7 1.1 5.6 0.6 1.7 0.5 1.2 1.0 0.2
## [46] 1.3 5.5 9.0 2.0 3.6data.frame(1:length(random.number),random.number,round(x,1))## X1.length.random.number. random.number round.x..1.
## 1 1 0.78809107 4.7
## 2 2 0.85113835 5.7
## 3 3 0.83712785 5.4
## 4 4 0.39579685 1.5
## 5 5 0.49387040 2.0
## 6 6 0.58318739 2.6
## 7 7 0.89667250 6.8
## 8 8 0.27145359 1.0
## 9 9 0.07705779 0.2
## 10 10 0.45359019 1.8
## 11 11 0.81436077 5.1
## 12 12 0.17863398 0.6
## 13 13 0.65323993 3.2
## 14 14 0.10332750 0.3
## 15 15 0.78108581 4.6
## 16 16 0.63047285 3.0
## 17 17 0.88616462 6.5
## 18 18 0.94045534 8.5
## 19 19 0.15061296 0.5
## 20 20 0.77057793 4.4
## 21 21 0.29947461 1.1
## 22 22 0.95971979 9.6
## 23 23 0.25744308 0.9
## 24 24 0.63572680 3.0
## 25 25 0.55166375 2.4
## 26 26 0.90367776 7.0
## 27 27 0.49211909 2.0
## 28 28 0.02802102 0.1
## 29 29 0.90893170 7.2
## 30 30 0.15761821 0.5
## 31 31 0.99124343 14.2
## 32 32 0.25043783 0.9
## 33 33 0.41506130 1.6
## 34 34 0.60070053 2.8
## 35 35 0.44833625 1.8
## 36 36 0.14886165 0.5
## 37 37 0.21541156 0.7
## 38 38 0.31173380 1.1
## 39 39 0.84588441 5.6
## 40 40 0.17162872 0.6
## 41 41 0.43257443 1.7
## 42 42 0.15236427 0.5
## 43 43 0.32574431 1.2
## 44 44 0.28721541 1.0
## 45 45 0.07355517 0.2
## 46 46 0.34325744 1.3
## 47 47 0.83887916 5.5
## 48 48 0.95096322 9.0
## 49 49 0.48161121 2.0
## 50 50 0.69702277 3.6hist(x, breaks = 20, probability = TRUE, col = "lightblue", ylim=c(0,0.35), main = "DistribuciĆ³n Exponencial (beta = 3)", xlab = "Valor de X")
# AƱadimos la funciĆ³n de densidad teĆ³rica de la distribuciĆ³n exponencial
curve(dexp(x,rate=1/3), col = "red", lwd = 2, add = TRUE)
Ejemplo: DistribuciĆ³n Poisson.
El nĆŗmero de veces que se cae la wifi de la universidad sigue una distribuciĆ³n Poisson con una media de 4 veces por horas. Simule el comportamiento de la caida de la red de la universidad.
\[ P(X=x)=\frac{e^{-4}4^x}{x!}\quad \text{para}\quad x=0,1,2,3,\dots \]
# Definir los parƔmetros
lambda <- 2 # Media de la distribuciĆ³n de Poisson
valores <- 0:7 # Valores de x de 0 a 10
# Inicializar vectores vacĆos para almacenar las probabilidades y las acumuladas
probabilidades <- numeric(length(valores))
acumuladas <- numeric(length(valores))
# Ciclo for para calcular la probabilidad de cada valor de x
for (x in valores) {
probabilidades[x + 1] <- dpois(x, lambda = lambda) # Probabilidad puntual
acumuladas[x + 1] <- ppois(x, lambda = lambda) # Probabilidad acumulada
}
# Mostrar los resultados
for (x in valores) {
cat("P(X =", x, ") =", probabilidades[x + 1], "\t Acumulada P(X <=", x, ") =", acumuladas[x + 1], "\n")
}## P(X = 0 ) = 0.1353353 Acumulada P(X <= 0 ) = 0.1353353
## P(X = 1 ) = 0.2706706 Acumulada P(X <= 1 ) = 0.4060058
## P(X = 2 ) = 0.2706706 Acumulada P(X <= 2 ) = 0.6766764
## P(X = 3 ) = 0.180447 Acumulada P(X <= 3 ) = 0.8571235
## P(X = 4 ) = 0.09022352 Acumulada P(X <= 4 ) = 0.947347
## P(X = 5 ) = 0.03608941 Acumulada P(X <= 5 ) = 0.9834364
## P(X = 6 ) = 0.0120298 Acumulada P(X <= 6 ) = 0.9954662
## P(X = 7 ) = 0.003437087 Acumulada P(X <= 7 ) = 0.9989033La regla para generar esta variable aleatoria estarĆa definida por:
[ { \[\begin{aligned} 0 & \quad \text{si}\, u_i \leq 0.1353\\ 1 & \quad \text{si}\, u_i \in (0.1353 , 0.4060] \\ 2 & \quad \text{si}\, u_i \in (0.4060 , 67667]\\ 3 & \quad \text{si}\, u_i \in (0.6767 , 0.8571] \\ 4 & \quad \text{si}\, u_i \in (0.8571, 0.9473]\\ 5 & \quad \text{si}\, u_i \in (0.9473 , 0.9834] \\ 6 & \quad \text{si}\, u_i \in (0.9834,0.9954]\\ 7 & \quad \text{si}\, u_i \in (0.9954 , 0.9999] \end{aligned}\]. ]
Algoritmos predefinidos.
Muchas librerĆas de software estadĆstico y de programaciĆ³n como NumPy en Python, R, MATLAB, y otros tienen funciones predefinidas para generar variables aleatorias de diversas distribuciones (normal, uniforme, exponencial, binomial, etc.), facilitando el proceso. Algunas de estas funciones utilizan mĆ©todos como el de la transforma inversa , convoluciĆ³n , composiciĆ³n , etc, .
SimulaciĆ³n de variables aleatorias Uniformes en R.
# Generar 100 variables aleatorias uniformes entre 0 y 1
uniform_random_01 <- runif(100, min = 0, max = 1)
uniform_random_01## [1] 0.4642812093 0.0472968640 0.4924413876 0.9845189112 0.2339436009
## [6] 0.6548774925 0.4724445250 0.8829573253 0.2922004333 0.6409208407
## [11] 0.3485153099 0.3120751686 0.2439736470 0.0001186794 0.6008252508
## [16] 0.8876450148 0.0191672409 0.8639476420 0.3889065925 0.9627142120
## [21] 0.4211694405 0.4152164590 0.9821176734 0.6708385306 0.1770130780
## [26] 0.1537170457 0.5387231887 0.1963973432 0.1832738505 0.8292213134
## [31] 0.6980703275 0.1290739689 0.8512812431 0.0164890783 0.6717730574
## [36] 0.3370865860 0.0310140641 0.2760329708 0.2162050249 0.0866783008
## [41] 0.6889543547 0.0289700108 0.0828108916 0.3146365935 0.1148433078
## [46] 0.4648038545 0.5544977824 0.1040754043 0.7515659479 0.2789934461
## [51] 0.8730775453 0.6871296221 0.9666396170 0.6632698372 0.3073290368
## [56] 0.3058602589 0.3791437685 0.8355505567 0.1519900281 0.3576006338
## [61] 0.8303941041 0.4876441127 0.8778959725 0.5569040247 0.0164652874
## [66] 0.5803387451 0.0395192457 0.2822865248 0.3296307956 0.1969050195
## [71] 0.5637719324 0.6260585375 0.5293665186 0.4855503265 0.1271492215
## [76] 0.5066240276 0.1703138223 0.9846078434 0.6361430045 0.3926924886
## [81] 0.0281733798 0.4750734994 0.4469799360 0.1050471738 0.0976128387
## [86] 0.7698865694 0.2296060033 0.8550958196 0.2901761879 0.8677754127
## [91] 0.2400990506 0.0416217335 0.4696948908 0.6046261955 0.9707362670
## [96] 0.2742241605 0.4766447518 0.9329083147 0.3595241059 0.1422693408# Generar 100 variables uniformes [ a, b]
uniform_random_ab <- runif(1000, min = 10, max = 20)Generando valores de una variable aleatoria normal con media \(\mu\) y varianza \(\sigma^2\) .
El mĆ©todo de la transformada inversa no es adecuado para la generaciĆ³n de variables aleatorias con distribuciĆ³n normal debido a la complejidad de la funciĆ³n de distribuciĆ³n acumulada (FDA) de la normal, cuya inversa no tiene una forma analĆtica simple. Esto hace que la aplicaciĆ³n directa de este mĆ©todo resulte impracticable en el caso de la distribuciĆ³n normal. Sin embargo, una alternativa eficaz es el uso del Teorema del LĆmite Central (TLC), el cual establece que la suma de un nĆŗmero suficientemente grande de variables aleatorias independientes e idĆ©nticamente distribuidas (i.i.d.) tiende a seguir una distribuciĆ³n normal, independientemente de la distribuciĆ³n original de estas variables, siempre que cumplan ciertas condiciones. Este enfoque permite aproximar la distribuciĆ³n normal mediante la suma de muchas variables mĆ”s simples, lo cual es una tĆ©cnica comĆŗn en simulaciones y mĆ©todos numĆ©ricos.
Teorema Central del LĆmite
Sea \(X_1, X_2, \dots, X_n\) una sucesiĆ³n de variables aleatorias independientes e idĆ©nticamente distribuidas (i.i.d.) con media \(\mu\) y varianza finita \(\sigma^2\). Entonces, la suma \(Y=\sum_{i=1}^{n} X_i\) tiende a seguir una distribuciĆ³n normal conforme \(n\) tiende a infinito. Es decir,
\[ Z_n = \frac{Y - n\mu}{\sigma \sqrt{n}} \sim N(0,1) \quad \text{cuando} \quad n \to \infty, \]
donde \(N(0,1)\) denota una distribuciĆ³n normal estĆ”ndar con media 0 y varianza 1.
En el caso de variables aleatorias uniforme \(R_i \sim iid U(0,1)\), entonces.
\[\begin{equation} Z = \frac{\sum_{i=1}^{n} R_i - \frac{n}{2}}{ \sqrt{\frac{n}{12}}} \sim N(0,1) \quad \text{cuando} \quad n \to \infty \label{\#eq:binom} \end{equation}\]
Puesto que la distribuciĆ³n normal estĆ”ndar de una variable aleatoria Y que se distribuye normal es \[Z=\frac{Y-\mu}{\sigma}\]
Por lo tanto la simulaciĆ³n de variables aleatorias \(Y\) se harĆa de acuerdo a la siguiente expresiĆ³n.
\[ Y=\mu +\sigma\times Z \]
con \(Z = \frac{\sum_{i=1}^{n} R_i - \frac{n}{2}}{ \sqrt{\frac{n}{12}}}\) Se ha comprobado que con \(n=12\) , la confiabilidad de los valores simulado es bastante aceptable, por lo que finalmente la expresiĆ³n nos queda.
\[ Y=\mu +\sigma\times \Big(\sum_{i=1}^{n} R_i -6\Big) \] Ejemplo:
# NĆŗmero de variables uniformes (cuanto mĆ”s grande n, mejor la aproximaciĆ³n)
n <- 12
# Generar 1000 variables aleatorias normales a travĆ©s de la convoluciĆ³n de uniformes
uniform_sums <- rowSums(matrix(runif(100 * n, min = 0, max = 1), nrow = 100, ncol = n))
uniform_sums## [1] 6.913273 4.091749 6.564155 7.652205 7.861741 6.841243 5.738543 5.287450
## [9] 7.739724 7.993367 6.648913 4.606619 7.880648 4.999721 6.533702 5.929031
## [17] 4.522285 5.882288 5.907364 7.354440 7.470344 5.302163 5.069417 8.084287
## [25] 6.906998 7.397485 6.829242 3.712715 7.000655 5.458774 5.970144 5.465639
## [33] 5.180264 4.971549 4.390170 6.049689 6.656151 5.221446 6.204418 5.374643
## [41] 4.891508 6.620499 6.315845 6.722656 7.603559 6.849221 6.962748 5.400733
## [49] 5.703208 6.215080 6.881678 4.549412 4.008492 5.278005 4.501137 6.598123
## [57] 5.120814 4.060961 6.658052 5.445219 5.512440 5.883089 6.334766 7.380367
## [65] 6.166718 7.433626 6.912407 5.107523 6.066910 4.873801 5.084337 7.090846
## [73] 6.297493 7.490770 4.508976 7.357550 5.240373 5.791605 5.400084 6.573421
## [81] 4.677110 5.129990 5.017883 5.954246 5.528616 6.742700 6.654290 7.630049
## [89] 5.097980 3.990192 6.360197 5.732873 6.941665 7.026370 4.576328 4.250700
## [97] 5.902974 5.514491 3.532195 3.892253# NormalizaciĆ³n para obtener una distribuciĆ³n aproximadamente N(0, 1)
Z <- (uniform_sums - 6)
# Mostrar los resultados
hist(Z, main = "AproximaciĆ³n de distribuciĆ³n normal usando convoluciĆ³n",
xlab = "Valores", col = "lightblue", probability = TRUE)
lines(density(Z), col = 2, lwd = 2)
# Simulando el tiempo de vida de 100 baterias de una lapto cuyo promedio de vida es de 30 meses y desviaciĆ³n tĆpica $\sigma=3.5$ meses, suponiendo que sigue una distribuciĆ³n normal
mu <- 30
sigma <- 3.5
Y= mu +3.5*Z
round(Y,1) # redondeando a una cifra decimal## [1] 33.2 23.3 32.0 35.8 36.5 32.9 29.1 27.5 36.1 37.0 32.3 25.1 36.6 26.5 31.9
## [16] 29.8 24.8 29.6 29.7 34.7 35.1 27.6 26.7 37.3 33.2 34.9 32.9 22.0 33.5 28.1
## [31] 29.9 28.1 27.1 26.4 24.4 30.2 32.3 27.3 30.7 27.8 26.1 32.2 31.1 32.5 35.6
## [46] 33.0 33.4 27.9 29.0 30.8 33.1 24.9 23.0 27.5 24.8 32.1 26.9 23.2 32.3 28.1
## [61] 28.3 29.6 31.2 34.8 30.6 35.0 33.2 26.9 30.2 26.1 26.8 33.8 31.0 35.2 24.8
## [76] 34.8 27.3 29.3 27.9 32.0 25.4 27.0 26.6 29.8 28.4 32.6 32.3 35.7 26.8 23.0
## [91] 31.3 29.1 33.3 33.6 25.0 23.9 29.7 28.3 21.4 22.6hist(Y, main = "Vida media de 100 baterias de laptos", xlab = "Valores", col = "lightblue", probability = TRUE)
#yi <- which(Y<24)
#yi
yi <- ifelse(Y<24,1,0)
prob <- sum(yi)/length(Y) # calculando la probabilidad impirica de que la vida media de una lapto sea inferior a 2 aƱos
prob## [1] 0.08Con el programa R funciĆ³n rnorm()
# Generar 100 variables normales con media 30 y desviaciĆ³n estĆ”ndar 3.5
normal_random <- rnorm(100, mean = 30, sd = 3.5)
normal_random## [1] 31.74991 26.03443 32.20652 27.83095 33.45688 33.15582 31.90606 32.09085
## [9] 32.81554 31.17961 24.48350 32.62965 27.35459 28.38254 29.16989 26.79850
## [17] 32.55509 36.78731 35.86258 31.55637 28.61217 29.79844 27.00242 27.87620
## [25] 31.03213 35.26591 30.37368 28.14310 27.31493 30.27267 24.59014 26.44837
## [33] 30.03863 28.23269 30.55801 31.06121 25.54045 19.97464 29.68502 32.69777
## [41] 24.87553 30.23975 32.42830 27.64673 32.35313 32.51736 33.82468 26.76639
## [49] 29.13182 27.48387 30.19282 28.57452 36.09113 28.28638 29.37664 27.40552
## [57] 29.66222 33.33518 30.68307 30.65393 31.10008 30.87994 34.68114 30.76478
## [65] 30.06513 35.77008 33.51283 24.32172 27.03253 22.82419 35.02900 35.13150
## [73] 29.19236 28.22154 29.61559 32.65524 35.46680 26.23591 31.50208 29.28553
## [81] 24.70775 26.79476 32.83444 30.93522 28.11897 27.80038 33.21022 26.90146
## [89] 30.99628 30.43192 30.18413 30.41068 30.70655 35.69003 27.29126 32.68519
## [97] 26.02514 26.42199 33.40020 25.85354media <- mean(Y)
v <- length(Y)-1
t <- qt(0.025,v,lower.tail=FALSE)
s <- sd(Y)
ee <- t*s/sqrt(length(Y))
lim_inf <- media-ee
lim_sup <- media + ee
cbind(lim_inf,lim_sup)## lim_inf lim_sup
## [1,] 29.03923 30.5906## [1] 0.02DistribuciĆ³n Exponencial
Para generar variables aleatorias con distribuciĆ³n exponencial con
parƔmetro \(\beta\).
Ejemplo:
Simular El tiempo de servicio en una caja de un banco cuyo comportamiento sigue una forma exponencial con media de 3 minutos/cliente
numexp <- rexp(100,1/3)
round(numexp,1)## [1] 0.4 2.2 1.6 0.9 0.6 9.0 1.1 0.7 5.9 0.8 1.1 0.2 0.1 2.6 1.6
## [16] 0.4 5.6 0.0 0.1 5.7 0.3 3.1 0.2 0.4 0.0 2.0 0.4 11.4 7.9 1.2
## [31] 5.2 0.8 3.2 2.9 3.8 0.2 9.3 0.8 0.5 1.1 5.7 2.2 1.3 4.3 1.3
## [46] 0.9 0.9 0.3 8.3 1.8 1.0 0.3 1.7 1.6 2.8 0.5 0.5 0.4 5.7 2.9
## [61] 0.5 6.9 6.0 0.7 4.3 6.8 0.4 1.3 2.3 0.2 4.9 2.3 8.9 0.2 0.3
## [76] 7.2 5.5 2.3 1.4 1.9 3.5 0.4 0.1 1.2 2.5 3.5 3.5 1.2 1.1 6.9
## [91] 1.0 8.7 2.4 0.6 0.5 3.6 2.8 2.0 2.2 0.9hist(numexp, breaks = 10, probability = TRUE, col = "lightblue", ylim=c(0,0.35), main = "DistribuciĆ³n Exponencial (beta = 3)", xlab = "Valor de X")
curve(dexp(x,rate=1/3), col = "red", lwd =2,add = TRUE)
DistribuciĆ³n Binomial
Simulando una prueba donde se realiza 12 preguntas con 4 opciones y una Ćŗnica respuesta, los estudiante responden de forma aleatoria.
prueba <- rbinom(1000000,12,0.25)
tabla_frec <- table(prueba)/length(prueba)
tabla_frec## prueba
## 0 1 2 3 4 5 6 7
## 0.031051 0.127398 0.232061 0.257740 0.194333 0.102958 0.040212 0.011464
## 8 9 10 11
## 0.002365 0.000378 0.000039 0.000001ganan <- ifelse(prueba>6,1,0)
sum(ganan)/length(prueba)## [1] 0.014247barplot(tabla_frec,col=3)
sum(tabla_frec)## [1] 1opciones <- c("A","B","C","D")
preguntas <- 1:12
respuestas <- sample(opciones,12,replace=TRUE)
cbind(preguntas,respuestas)## preguntas respuestas
## [1,] "1" "C"
## [2,] "2" "D"
## [3,] "3" "B"
## [4,] "4" "D"
## [5,] "5" "D"
## [6,] "6" "B"
## [7,] "7" "D"
## [8,] "8" "C"
## [9,] "9" "B"
## [10,] "10" "B"
## [11,] "11" "B"
## [12,] "12" "B"DistribuciĆ³n Poisson
Supongamos que queremos simular un proceso donde los eventos ocurren con una tasa de \(\lambda=5\) eventos por minuto .
evento_pois <- rpois(100,5)
tablapois <- table(evento_pois)
tablapois## evento_pois
## 1 2 3 4 5 6 7 8 9 10
## 2 7 16 21 19 17 6 5 4 3# GrƔfico de los eventos
barplot(tablapois,col = "blue")
DistribuciĆ³n Erlang
La distribuciĆ³n Erlang es un modelo ampliamente utilizado en sistemas que involucran tiempos de espera y procesos en fases, como en la teorĆa de colas, redes de telecomunicaciones y anĆ”lisis de supervivencia. La distribuciĆ³n Erlang es Ćŗtil para modelar el tiempo hasta que ocurren k eventos en un proceso de Poisson, donde cada evento tiene un tiempo de espera exponencial independiente con el mismo parĆ”metro de tasa \(\lambda\), es un caso especial de la distribuciĆ³n Gamma cuando el parĆ”metro de forma \(k\) es entero
\[f(x; k, \lambda) = \frac{\lambda^k x^{k-1} e^{-\lambda x}}{(k-1)!}, \quad x \geq 0\]
# simulando variables aleatoria Erlang : "Tres llamadas a un tasa de dos minutos"
k=3
lambda=1
Erl<-rgamma(100,shape=k,rate=lambda)
hist(Erl,col=6,probability=TRUE)
lines(density(Erl),col="blue",lty=1)
Erl## [1] 0.4666872 1.3733782 3.8021897 2.8666580 2.2115104 4.9890526 1.9321328
## [8] 1.9662018 2.3376292 0.6879687 2.0533757 3.2187923 2.1752124 4.2833493
## [15] 1.3910484 2.2737018 0.6623244 2.4790170 2.8267075 1.3620061 5.0202938
## [22] 0.4508373 0.9797896 3.8664258 2.0084628 2.3169570 1.8520915 3.6466546
## [29] 1.7049884 1.1205509 5.0806210 4.7918321 4.3466045 2.9264931 2.9299101
## [36] 0.8996342 5.7169492 2.1636658 4.4391427 2.4914176 1.1754977 1.4381781
## [43] 2.6586577 1.5231518 2.0248956 3.1276793 1.9169547 1.5631375 3.3078683
## [50] 3.6465393 0.9607019 6.4805216 1.7844781 4.6231707 5.9391728 3.0228432
## [57] 2.8657089 1.4464842 1.0244442 4.6562605 3.5224538 4.4789741 3.6363859
## [64] 5.5446336 2.4155124 2.5182384 0.7540953 1.9542661 3.2785437 3.9121243
## [71] 1.8472351 5.2004369 0.7391280 3.4935333 3.5174236 4.9642440 2.5347602
## [78] 6.5399099 2.2901547 5.1292132 2.9082922 1.8525222 2.0609848 2.8783444
## [85] 6.1331680 0.8302160 3.7694358 2.6691781 1.1092151 4.1867868 3.9118532
## [92] 3.5417205 0.7616364 3.5107945 0.7415535 2.9239951 2.9747245 2.2314757
## [99] 3.3774329 3.6697998Ejemplo de Aplicaciones de SimulaciĆ³n.
Modelos Terminales
Los modelos terminales de simulaciĆ³n se caracterizan por un evento que pone fin a la simulaciĆ³n. Un ejemplo serĆa medir el tiempo necesario para procesar 10 piezas, vender 10000 pc o atender clientes en una cafeterĆa en un perĆodo determinado. El anĆ”lisis estadĆstico recomendado incluye el uso de intervalos de confianza y la identificaciĆ³n de la distribuciĆ³n de probabilidad de la variable de salida.
Modelos no terminales Las simulaciones no terminales o de estado estable no tienen un punto de finalizaciĆ³n definido en el tiempo. Por ejemplo, al modelar la cantidad de mĆ”quinas necesarias en un sistema de producciĆ³n que opera continuamente, la simulaciĆ³n continĆŗa hasta que la variable de interĆ©s alcance un estado estable. Es importante determinar la longitud de la corrida para garantizar que los resultados del modelo se estabilicen adecuadamente.
Modelo de una lĆnea de espera con un servidor
El tiempo que transcurre entre la llegada de ciertas piezas a una estaciĆ³n de inspecciĆ³n sigue una distribuciĆ³n exponencial con media de 5 minutos/pieza. El proceso estĆ” a cargo de un operario, y la duraciĆ³n de la inspecciĆ³n sigue una distribuciĆ³n normal con media de 6.0 y desviaciĆ³n estĆ”ndar de 1.0 minutos/pieza. Calcule el tiempo promedio de permanencia de las piezas en el proceso de inspecciĆ³n.
piezas <- 1:1000 # NĆŗmero de piezas simuladas
media_llegadas <- 5 # parƔmetro de la exponencial (para el tiempo entre llegadas)
media_inspeccion <- 6 # Media para el tiempo de inspecciĆ³n (distribuciĆ³n normal)
sd_inspeccion <- 1 # DesviaciĆ³n estĆ”ndar
# Generar tiempos entre llegadas (distribuciĆ³n exponencial)
tiempos_llegadas<-rexp(n,rate=1/media_llegadas)
# Generar tiempos de inspecciĆ³n (distribuciĆ³n normal)
tiempos_inspeccion <- rnorm(n, mean = media_inspeccion, sd = sd_inspeccion)
# Inicializar el tiempo de llegada y salida de cada pieza
tiempo_actual <- 0
tiempos_perm <- numeric(n) # Para almacenar los tiempos de permanencia
for (i in 1:n) {
tiempo_llegada <- tiempo_actual + tiempos_llegadas[i]
tiempo_salida <- tiempo_llegada + tiempos_inspeccion[i]
# Calcular el tiempo de permanencia en el sistema
tiempos_perm[i] <- tiempo_salida - tiempo_actual
# Actualizar el tiempo actual al tiempo de salida
tiempo_actual <- tiempo_salida
}
# Calcular el tiempo promedio de permanencia
tiempo_promedio_perm <- mean(tiempos_perm)
# Mostrar resultado
cat("El tiempo promedio de permanencia en el proceso de inspecciĆ³n es:", tiempo_promedio_perm, "minutos.\n")## El tiempo promedio de permanencia en el proceso de inspecciĆ³n es: 10.81556 minutos.