Simulación de distribuciones continuas

Para detalles de las distribuciones puede obtener el libro “Modelos de probabilidad” en la siguiente dirección

https://1drv.ms/b/s!Aj-hHTVbsx01hs9cx6tbr9v1vZKjzQ?e=qyn6U4

1. Método de la transformada inversa.

El método de la transformada inversa se utiliza para simular variables aleatorias continuas, mediante la función de distribución acumulada y la generación de números aleatorios con distribución uniforme en el intervalo (0,1).

1.1. Distribución uniforme.

La función de densidad de probabilidad de una variable aleatoria continua con distribución uniforme en un intervalo [a,b]:

\[\begin{equation*} f(x)= \begin{cases} \dfrac{1}{b-a} &\text{si $a \leq x \leq b$}\\ 0 &\text{en otro caso} \end{cases} \end{equation*}\]

La función de distribución acumulada de la variable aleatoria continua con distribución uniforme en el intervalo

\[\begin{equation*} F(x)= \begin{cases} 0 &\text{si $x < a$}\\ \dfrac{x-a}{b-a} &\text{si $a \leq x \leq b$} \end{cases} \end{equation*}\]

Se hace: \(u=\dfrac{x-a}{b-a}\)

\(x-a=u(b-a)\) \[\begin{equation*} x=a+(b-a)u \end{equation*}\]

Ejemplo. Simular una variable uniforme continua en el intervalo [50,70]

a<-50;b<-70
u<-runif(100000,0,1)
x<-a+(b-a)*u
hist(x,prob=TRUE,main="",ylab="f(x)",col="grey",ylim=c(0,0.06))
curve(dunif(x,a,b),  add = TRUE,lwd=2)
box()
text(60,0.055,expression(paste(a==50, b==70)))

1.2. Distribución exponencial.

Sea una variable aleatoria continua, $ X$, que tiene función de densidad: \[\begin{equation*} f(x)= \begin{cases} \dfrac{1}{\theta} e^{-\frac{x}{\theta}} &\text{si $ x \geq 0$}\\ 0 &\text{en otro caso} \end{cases} \end{equation*}\]

La función de distribución acumulada de la distribución exponencial es: \[\begin{equation*} F(x)= \begin{cases} 0 & \text{si $ x < 0 $}\\ 1- e^{- \frac{x}{\theta} } &\text{si $ x \geq 0$}\\ \end{cases} \end{equation*}\]

Se hace:

\(u=1-e^{- \frac{x}{\theta} }\)

\(1-u=e^{- \frac{x}{\theta} }\)

\(ln\,(1-u)= - \frac{x}{\theta}\)

\[\begin{equation*} x= - \theta \cdot ln(1-u) \end{equation*}\]

Ejemplo. Simular una variable aleatoria con distribución exponencial con \(\theta =2\).

theta<-2
u<-runif(100000,0,1)
x<--theta*log(1-u)
hist(x,prob=TRUE,ylab="f(x)",main="",col="grey")
curve(1/theta*exp(-x/theta),add=TRUE,lwd=2)
box()
text(10,0.25,expression(paste(theta==2.0)))

1.3. Distribución Weibull.

La función de densidad de probabilidad de una variable aleatoria con distribución Weibull con parámetros \(\alpha\) y \(\beta\) es:

\[\begin{equation*} f(x)= \begin{cases} \dfrac{\alpha}{\beta} \left(\dfrac{x}{\beta} \right)^{\alpha - 1}e^{-\left( \dfrac{x}{\beta} \right)^{\alpha} }&\text{si $ x >0$}\\ 0 &\text{en otro caso} \end{cases} \end{equation*}\]

La función de distribución acumulada de la distribución Weibull es:

\[\begin{equation*} F(x)= \begin{cases} 0 &\text{si $ x <0$}\\ 1-e^{- \left( \dfrac{x}{\beta} \right)^{\alpha} } &\text{si $ x \geq 0$} \end{cases} \end{equation*}\]

Si se hace:

\(u= 1-e^{- \left( \dfrac{x}{\beta} \right)^{\alpha} }\)

\(1-u=e^{- \left( \dfrac{x}{\beta} \right)^{\alpha} }\)

\(ln\,(1-u)= - \left(\dfrac{x}{\beta} \right)^{\alpha}\)

\(x^{\alpha}= - \beta^{\alpha}\cdot ln(1-u)\)

\[\begin{equation*} x= \sqrt[\alpha]{ - \beta ^{\alpha} \cdot ln\,(1-u)} \end{equation*}\]

Ejemplo. Simular una variable aleatoria que tiene distribución Weibull con parámetros \(\alpha= 6\) y \(\beta= 0.5\)

alpha<-6;beta<- 0.5
u<-runif(100000,0,1)
x=(-beta^alpha*log(1-u))^(1/alpha)
hist(x,prob=TRUE,ylab="f(x)", main="",ylim=c(0,5),col="grey")
curve(dweibull(x,alpha,beta),add=TRUE,lwd=2)
box()
text(0.2,3,expression(paste(alpha==6)))
text(0.2,2.7,expression(paste(beta==0.5)))

1.4. Distribución Gumbel.

La función de densidad de probabilidad de una variable aleatoria con distribución Gumbel con parámetros \(\mu\) y \(\beta\) es: \[\begin{equation*} f(x)=\dfrac{1}{\beta}exp\left[ - \left( \dfrac{x-\mu}{\beta} \right) exp\left( - \left( \dfrac{x-\mu}{\beta} \right) \right)\right] \end{equation*}\]

La función de distribución acumulada es:

\[\begin{equation*} F(x)=exp\left( -exp\left(-\dfrac{x-\mu}{\beta} \right) \right) \end{equation*}\]

Si se hace:

\[\begin{equation*} x= \mu-\beta\,ln\left( ln\,\left(\dfrac{1}{u} \right) \right) \end{equation*}\]

Ejemplo. Simular una distribución Gumbel con parámetros \(\mu= 20\) y \(\beta= 5\).

u<-runif(100000,0,1)
mu<-20
beta<-5
x<-mu-beta*log(log(1/u))
hist(x,prob=T,breaks="Sturges",ylab="f(x)",xlab="x",main="",col="grey")
par(new=T)
y<-function(x){
1/beta*exp(-(x-mu)/beta)*exp(-exp(-(x-mu)/beta))
}
curve(y,10,50,lwd=2,axes=F,ylab="")
box()
text(40,0.05,expression(paste(mu==20)))
text(40,0.045,expression(paste(beta==5)))

1.5. Distribución Fréchet.

\[\begin{equation*} f(x)=\dfrac{\alpha}{\delta}\left( \dfrac{x-\lambda}{\delta} \right)^{-1 - \alpha} \cdot exp\left( -\left( \dfrac{x - \lambda}{\delta}\right)^{- \alpha} \right) \end{equation*}\]

Para \(x >\lambda\)

\(\alpha\): parámetro de forma, \(\delta\): parámetro de escala y \(\lambda\): parámetro de localización.

\[\begin{equation*} F(x)= exp\left( - \left( \dfrac{x- \lambda}{\delta}\right)^{- \alpha} \right) \end{equation*}\] Para \(x > \lambda\)

Se hace

\(u=exp\left( - \left( \dfrac{x- \lambda}{\delta}\right)^{- \alpha} \right)\)

\(ln\,(u)= -\left( \dfrac{x-\lambda}{\delta} \right)^{-\alpha}\)

\(\left[ ln\,\left( \dfrac{1}{u}\right)\right] ^{-\dfrac{1}{\alpha}} =\dfrac{x-\lambda}{\delta}\)

\(x=\lambda+\delta\left[ ln\,\left( \dfrac{1}{u} \right) \right]^{- \dfrac{1}{\alpha}}\)

Ejemplo. Simular una distribución Fréchet con parámetros \(\alpha= 7\), \(\delta=4\) y \(\lambda= 20\).

u<-runif(10000,0,1)
l<-20;d<-4;a<-7
x<-l+d*(-log(u))^(-1/a)
hist(x,prob=T,col="grey",main="",ylab="f(x)")
y<-function(x){
(a/d)*((x-l)/d)^(-1-a)*exp((-(x-l)/d)^(-a))
}
par(new=T)
curve(y,min(x),30,axes=FALSE,lwd=2,ylab="")
box()
text(27,0.43,expression(paste(lambda==20)))
text(27,0.40,expression(paste(delta==4)))
text(27,0.37,expression(paste(alpha==7)))

1.6. Distribución Pareto.

La función de densidad de una variable aleatoria continua con distribución Pareto tiene la forma: \[\begin{equation*} f(x)= \begin{cases} \dfrac{\alpha}{x}\left( \dfrac{\beta}{x}\right) ^{\alpha} &\text{si $ x \geq \beta$}\\ $0$ &\text{en otro caso} \end{cases} \end{equation*}\] La función de distribución acumulada es: \[\begin{equation*} F(x)= \begin{cases} 0 &\text{si $ x < \beta$}\\ 1-\left( \dfrac{\beta}{x} \right)^{\alpha} &\text{si $ x \geq \beta$} \end{cases} \end{equation*}\]

Se hace:

\[\begin{equation*} x= exp\left( {\frac{\alpha\cdot ln \beta - ln (1- u)}{\alpha}}\right) \end{equation*}\]

Ejemplo. Sea una variable aleatoria con distribución Pareto de parámetros \(\alpha = 150\) y \(\beta = 25.5\)

beta=150
alpha=25.5
u<-runif(10000,0,1)
x<-exp((alpha*log(beta)-log(1-u))/alpha)
hist(x,prob=TRUE,main="",ylab="f(y)", ylim=c(0,0.10),col="grey")
y<-function(x){
(alpha/x)*(beta/x)^(alpha)
}
par(new=T)
curve(y,min(x),200,axes=FALSE,lwd=2,ylab="")
box()
text(180,0.1,expression(paste(beta==150)))
text(180,0.09,expression(paste(alpha==25.5)))

1.7. Distribución logística.

La función de densidad de una variable aleatoria continua con distribución Pareto tiene la forma: \[\begin{equation*} \dfrac{exp\left[- \left( \dfrac{x-\alpha}{\beta}\right) \right] }{\beta\left\lbrace1+exp\left[ -\left(\dfrac{x-\alpha}{\beta} \right) \right] \right\rbrace ^{2}} \end{equation*}\]

La función de distribución acumulada es: \[\begin{equation*} F(x)=\dfrac{1}{1+exp\left[ \left( -\dfrac{x-\alpha}{\beta}\right) \right] } \end{equation*}\]

Se hace:

\[\begin{equation*} x=\alpha- \beta \cdot ln\left(\dfrac{1-u}{u} \right) \end{equation*}\]

Ejemplo. Sea una variable aleatoria con distribución logística de parámetros \(\alpha = 100\) y \(\beta = 15\)

alpha=100;beta=15
u<-runif(100000,0,1)
x<-alpha-beta*log((1-u)/u)
hist(x,prob=T,main="",ylab="f(x)",ylim=c(0,0.016),col="grey")
y<-function(x){
exp(-(x-alpha)/beta)/(beta*( 1+exp(-(x-alpha)/beta))^2)
}
par(new=T)
curve(y,-50,250,main="",axes=FALSE,lwd=2,ylab="")
box()
text(180,0.01,expression(paste(alpha==100)))
text(180,0.009,expression(paste(beta==15)))

2. Método de convolución.

En algunas distribuciones de probabilidad la variable aleatoria a simular puede generarse mediante la suma de otras variables aleatorias, de manera tal que es más rápida la simulación que utilizando otros métodos.

El método de convolución se puede expresar de la siguiente forma: \[\begin{equation*} Y=X_{1}+X_{2}+X_{3}+\cdots+X_{k} \end{equation*}\]

2.1. Distribución gamma.

La función de densidad de una variable aleatoria continua con distribución gamma con parámetros \(\alpha\) y \(\beta\) es: \[\begin{equation*} f(x)= \begin{cases} \dfrac{1}{\Gamma(\alpha)\beta^{\alpha}}x^{\alpha-1}e^{-\dfrac{x}{\beta}} &\text{si $ x \geq 0$}\\ 0 &\text{en otro caso} \end{cases} \end{equation*}\]

Para simular una variable aleatoria con distribución gamma de parámetros \(\alpha\) y \(\beta\), donde \(\alpha\) es un entero, se utiliza el hecho de que la suma de \(n\) variables aleatorias independientes con distribución exponencial de parámetro $ $ es una variable aleatoria con distribución gamma.

Entonces, si \(U_{1}\), \(U_{2}\), \(U_{3}\), \(\ldots\), \(U_{n}\) son variables aleatorias independientes con distribución uniforme en \((0,1)\) \[\begin{equation*} x=-\displaystyle\sum_{i=1}^{n}\dfrac{1}{\theta}\,ln\,u_{i}= -\dfrac{1}{\theta} ln\,\left( \displaystyle\prod _{i=1}^{n}u_{i}\right) \end{equation*}\]

Ejemplo. Simular una variable aleatoria continua con distribución gamma de parámetros \(\alpha=4\) y \(beta=2.8\).

alpha=4;beta=2.8
u1=runif(100000,0,1);u2=runif(100000,0,1)
u3=runif(100000,0,1);u4=runif(100000,0,1)
x=-beta*log(u1*u2*u3*u4)
hist(x,prob=TRUE,breaks=9,ylim=c(0,0.09),main="" ,ylab="f(x)",col="grey")
y<-function(x){
1/(gamma(alpha)*beta^alpha)*x^(alpha-1)*exp(-x/beta)
}
par(new=T)
curve(y,0,50,main="",axes=FALSE,lwd=2,ylab="")
box()
text(30,0.07,expression(paste(beta==2.8)))
text(30,0.065,expression(paste(alpha==4)))

2.2. Distribución normal.

La variable aleatoria normal con media \(\mu\) y varianza \(\sigma^{2}\) se puede generar utilizando el .

\(Y=X_{1}+X_{2}+X_{3}+\cdots+X_{k}\) \(\sim N \left( k \mu_{X}, k \sigma_{X}^{2}\right)\).

Al reemplazar \(X\) por números pseudoaleatorios,

\(Y=U_{1}+U_{2}+U_{3}+\cdots +U_{k} \sim N \left( k\dfrac{1}{2},k\dfrac{1}{12}\right)\).

Ejemplo. Simular una variable aleatoria normal con media 3 y varianza 0.5.

u1=runif(100000,0,1);u2=runif(100000,0,1)
u3=runif(100000,0,1);u4=runif(100000,0,1)
u5=runif(100000,0,1);u6=runif(100000,0,1)
y=u1+u2+u3+u4+u5+u6
hist(y,prob=TRUE,breaks=15,ylim=c(0,0.6),main="", ylab="f(x)",col="grey")
curve(dnorm(x, mean=3,sd=0.7071), add = TRUE,lwd=2,ylim=c(0,0.6))
box()
text(5,0.5,expression(paste(mu==3)))
text(5,0.45,expression(paste(sigma==0.5)))

Ejemplo. Utilizar el método anterior para generar una variable aleatoria normal con media 6 y varianza 1.

u1=runif(100000,0,1);u2=runif(100000,0,1);u3=runif(100000,0,1)
u4=runif(100000,0,1);u5=runif(100000,0,1);u6=runif(100000,0,1)
u7=runif(100000,0,1);u8=runif(100000,0,1);u9=runif(100000,0,1)
u10=runif(100000,0,1);u11=runif(100000,0,1);u12=runif(100000,0,1)
y=u1+u2+u3+u4+u5+u6+u7+u8+u9+u10+u11+u12
hist(y,prob=TRUE,breaks=7,ylim=c(0,0.5),main="",ylab="f(x)",col="grey")
curve(dnorm(x, mean=6,sd=1), add = TRUE,lwd=2,ylim=c(0,0.5))
box()
text(9,0.4,expression(paste(mu==6)))
text(9,0.35,expression(paste(sigma==1)))

2.3. Normal estándar.

A partir de la anterior distribución se puede generar la distribución normal estándar, haciendo

\(Z= U_{1}+U_{2}+U_{3}+U_{4}+U_{5}+U_{6}+U_{7}+U_{8}+U_{9}+U10+U_{11}+U_{12}-6 \sim N\left( 0,1\right)\).

Ejemplo. Utilizar el procedimiento anterior para generar una variable normal estándar.

u1=runif(100000,0,1);u2=runif(100000,0,1);u3=runif(100000,0,1)
u4=runif(100000,0,1);u5=runif(100000,0,1);u6=runif(100000,0,1)
u7=runif(100000,0,1);u8=runif(100000,0,1);u9=runif(100000,0,1)
u10=runif(100000,0,1);u11=runif(100000,0,1);u12=runif(100000,0,1)
z=(u1+u2+u3+u4+u5+u6+u7+u8+u9+u10+u11+u12)-6
hist(z,prob=TRUE,breaks=15,ylim=c(0,0.4),main="",ylab="f(x)",col="grey")
curve(dnorm(x,mean=0,sd=1), add = TRUE,lwd=2,ylim=c(0,0.4))
box()
text(3,0.35,expression(paste(mu==0)))
text(3,0.32,expression(paste(sigma==1)))

2.4. Distribución normal.

De la simulación de la variable con distribución normal estándar:

\(Z= U_{1}+U_{2}+U_{3}+U_{4}+U_{5}+U_{6}+U_{7}+U_{8}+U_{9}+U10+U_{11}+U_{12}-6 \sim N\left( 0,1\right)\).

\(Z=\displaystyle\sum_{i=1}^{12}u_{i}-6\)

\(\dfrac{x-\mu}{\sigma}=\displaystyle\sum_{i=1}^{12}u_{i}-6\)

Despejando \(x\), se tiene que :

\[\begin{equation*} x=\left[ \displaystyle\sum_{i=1}^{n}u_{i}-6\right] \sigma+\mu \end{equation*}\] que permite generar una variable aleatoria con distribución normal con media \(\mu\) y desviación estándar \(\sigma\).

Ejemplo. Generar una variable aleatoria continua que tiene distribución normal con media 40 y desviación estándar 7.

Se tiene que:

\[\begin{equation*} x=\left[ \displaystyle\sum_{i=1}^{n}u_{i}-6\right]7 +40 \end{equation*}\]

mu=40;sigma=7
u1=runif(100000,0,1);u2=runif(100000,0,1)
u3=runif(100000,0,1);u4=runif(100000,0,1)
u5=runif(100000,0,1);u6=runif(100000,0,1)
u7=runif(100000,0,1);u8=runif(100000,0,1)
u9=runif(100000,0,1);u10=runif(100000,0,1)
u11=runif(100000,0,1);u12=runif(100000,0,1)
x=((u1+u2+u3+u4+u5+u6+u7+u8+u9+u10+u11+u12)-6)*sigma+mu
summary(x)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   11.75   35.24   40.04   40.02   44.78   66.12

sd(x)

## [1] 6.996632

hist(x,prob=TRUE,breaks=15,ylim=c(0,0.06),main="", ylab="f(x)",col="grey")
curve(dnorm(x, mean=mu,sd=sigma), add = TRUE,lwd=2,ylim=c(0,0.06))
box()
text(60,0.04,expression(paste(mu==40)))
text(60,0.037,expression(paste(sigma==7)))

3. Método de transformación directa.

Este método se utiliza para generar variables aleatoria con distribución normal.

Ejemplo. Generar una variable aleatoria continua que tiene distribución normal estándar.

3.1. Normal estándar.

mu<-0;sigma<-1
u1<-runif(100000,0,1)
u2<-runif(100000,0,2*pi)
x1<- (sqrt(-2*log(1-u1))*cos(2*pi*u2))*sigma+mu
hist(x1,prob=TRUE,ylim=c(0,0.55), ylab="f(x)",col="grey",main=" Distribución normal estándar")
curve(dnorm(x, mean=mu,sd=sigma), add = TRUE,lwd=2,ylim=c(0,0.55))

3.2. Distribución normal.

Ejemplo. Generar una variable aleatoria continua que tiene distribución normal con media 50 y desviación estándar 6.

mu<-50;sigma<-6
u1<-runif(100000,0,1)
u2<-runif(100000,0,2*pi)
x1<- (sqrt(-2*log(1-u1))*cos(2*pi*u2))*sigma+mu
hist(x1,prob=TRUE,ylim=c(0,0.075), ylab="f(x)",col="grey",main=" Distribución normal estándar")
curve(dnorm(x, mean=mu,sd=sigma), add = TRUE,lwd=2,ylim=c(0,0.55))
box()

4. Metodo de composición.

A partir de este método se pueden generar variables aleatorias cuando la función de densidad se expresa como una combinación de \(m\) distribuciones, es decir, se encuentra definida por partes.

Para esta situación, la función de densidad se puede expresar como: \(f(x)=\displaystyle\sum_{i=1}^{m}f_{i}(x)I_{A}(x)\)

Donde \(I{A}(x)\) es la función indicadora definida como: \[\begin{equation*} I_{A}(x)= \begin{cases} 0 & \text{si $ x \in A $}\\ 1 & \text{si $x \notin A $} \end{cases} \end{equation*}\]

4.1. Distribución triangular.

La función de densidad de la variable aleatoria con distribución triangular es: \[\begin{equation*} f(x)= \begin{cases} \dfrac{2(x-a)}{(b-a)(c-a)} & \text{si $ a<x\leq c $}\\ \dfrac{2(b-x)}{(b-a)(b-c)} & \text{si $c<x\leq b $} \end{cases} \end{equation*}\] A continuación se calcula la probabilidad para cada una de las dos partes en que está definida la distribución. \[\begin{equation*} p(x)= \begin{cases} \displaystyle\int_{a}^{c}\dfrac{2(x-a)}{(b-a)(c-a)} dx =\dfrac{(c-a)}{(b-a)} & \text{si $a<x\leq c $}\\ \displaystyle\int_{c}^{b}\dfrac{2(b-x)}{(b-a)(b-c)} dx=\dfrac{(b-c}{(b-a)} & \text{si $c<x\leq b $} \end{cases} \end{equation*}\] Como cada uno de las partes separadamente no constituyen una función de probabilidad, se ajustan dividiendo por su correspondiente probabilidad. \[\begin{equation*} f(x)= \begin{cases} \dfrac{2(x-a)}{(b-a)(c-a)}\dfrac{(b-a)}{c-a)}=\dfrac{2(x-a)}{(c-a)^{2}} & \text{si $ a<x\leq c $}\\ \dfrac{2(b-x)}{(b-a)(b-c)}\dfrac{(b-a)}{(b-c)}=\dfrac{2(b-x)}{(b-c)^{2}} & \text{si $c<x\leq b $} \end{cases} \end{equation*}\] Para aplicar el método de la transformación inversa, se obtiene la función de distribución acumulada \[\begin{equation*} F(x)= \begin{cases} \displaystyle\int_{a}^{x} \dfrac{2(x-a)}{(c-a)^{2}}dx=\dfrac{(x-a)^{2}}{(c-a)^{2}} & \text{si $ a<x\leq c $}\\ \displaystyle\int_{c}^{x}\dfrac{2(b-x)}{(b-c)^{2}}dx=1-\dfrac{(b-x)^{2}}{(b-c)^{2}} & \text{si $c<x\leq b $} \end{cases} \end{equation*}\] Despejando \(x\) y reemplazando por los números pseudoaleatorios, en la función de distribución acumulada, se obtiene: \[\begin{equation*} x= \begin{cases} a+(c-a)\sqrt{u_{i}} & \text{si $ a<x\leq c $}\\ b-(b-c)\sqrt{1-u_{i}} & \text{si $c<x\leq b $} \end{cases} \end{equation*}\] Expresando la anterior ecuación con la función indicadora, se obtiene: \[\begin{equation*} x= \begin{cases} a+(c-a)\sqrt{u_{i}} & \text{si $ u_{j}\leq \dfrac{(c-a)}{(b-a)} $}\\ b-(b-c)\sqrt{1-u_{i}} & \text{si $u_{j}>\dfrac{(c-a)}{(b-a)} $} \end{cases} \end{equation*}\]

Ejemplo. Generar una variable aleatoria continua con distribución triangular con parámetros \(a=10\), \(b=30\) y \(c=25\).

a<-10;b<-30;c<-25
ri<-(c-a)/(b-a)
u1<-runif(10000,0,1)
u2<-runif(10000,0,1)
for(x1 in u1 ){
if(x1<=ri){
x2<-a+((c-a)*sqrt(u1))
}else
x3<-b-((b-c)*sqrt(1-u2))
}
x=c(x2,x3)
hist(x,prob=TRUE,breaks=8, main="", ylab="f(x)",col="grey")
curve(dtriang(x, min=10, mode=25, max=30), from = 10, to = 30,lwd=2,add=TRUE) 
box()
text(16,0.1,expression(paste(a==40)))
text(16,0.09,expression(paste(b==30)))
text(16,0.08,expression(paste(c==25)))

5. Método de aceptación y rechazo.

Este método se genera una variable aleatoria continua \(Y\), que tiene una función de densidad conocida \(g(y)\), posteriormente se genera un número aleatorio \(U_{1}\) y se decide aceptar el valor generado con una probabilidad proporcional a \(\dfrac{f(x)}{g(y)}\). La función \(g(y)\) debe ser una función de de densidad que acote a \(f(x)\), pero en general, \(g(y)\) no es una función de densidad.

5.1. Distribución beta.

Ejemplo. Simular una variable aleatoria continua con distribución beta de parámetros \(\alpha=3\) y \(\beta=2\)

La función de densidad es:

\(f(x)=12x^{2}(1-x)\)

Esta función tiene un máximo en \(x=\dfrac{2}{3}\) y \(f(\frac{2}{3})=\dfrac{16}{9}\).

Se define \(g(x)=\dfrac{16}{9}\) para \(0\leq x \leq 1\). Entonces, \(g(x)\) es mayor que \(f(x)\) en (0,1)

Generar \(U_{1}\) y \(U_{2}\).

Si \(U_{2}\leq \dfrac{12 x^{2}(1-x)}{\frac{16}{9}}\)

hacer \(x=u_{1}\).

u1=runif(100000,0,1)
u2=runif(100000,0,1)
T=(12*u1^2*(1-u1))/(16/9)
x=u1[u2<=T]
n=length(x)
hist(x,prob=TRUE,breaks=15,ylim=c(0,2.1),main=" ",ylab="f(x)",col="grey")
curve(dbeta(x,3,2),lwd=2,add=TRUE)
box()
text(0.2,1.6,expression(paste(alpha==3)))
text(0.2,1.5,expression(paste(beta==2)))

|—|