Considere una variable aleatoria discreta \(X\) con función masa de probabilidad dada por \[ f_{X}(x)=\left(\begin{array}{c} n \\ x \end{array}\right) \frac{B(x+\alpha, n-x+\beta)}{B(\alpha, \beta)}, x=0,1,2, \ldots, n \] con \(\alpha>0, \beta>0\), \(n\) un número entero mayor que cero y \(B(\cdot, \cdot)\) la función beta. Esta distribución de probabilidad es conocida como distribución beta-binomial y corresponde a una distribución binomial en la que la probabilidad de éxito en cada uno de los \(n\) ensayos no es fija.
Cree una función que emplee el método de aceptación y rechazo, y use una distribución binomial como distribución de referencia, para generar una muestra de tamaño \(m= 10000\) de esta distribución (con \(n=25\), \(\alpha=2\) y \(\beta=5\) ) y grafique además su histograma junto con la función masa de probabilidad teórica (escriba en el pdf la función creada y la gráfica pedida).
# distribucion masa de probabilidad beta binomial
betabinom = function(x, size, A, B) {
choose(size,x) * beta(x+A, n-x+B) / beta(A,B)
}
# funcion generadora de numeros aleatorios beta binomial
g.betabinom = function(m, n, A, B){ # m: tamano de muestra
x = 0:n
c. = max(betabinom(x, size = n, A, B)/dbinom(x, size=n, prob = 0.432))
x.m = numeric(0) ; i = 0
while(length(x.m) < m){
j = rbinom(1, size = n, prob = 0.432) # paso 1
u = runif(1) # paso 2
if(u <= betabinom(j, size = n, A, B)/(c.*dbinom(j, size = n, prob = 0.432))){
i = i + 1
x.m[i] = j # paso 3
}
}
x.m
}
# Parametros para la grafica
n = 25
A = 2
B = 5
X = g.betabinom(m = 1000
, n, A, B)
plot(table(X)/sum(table(X)))
lines(0:n, betabinom(0:n, n, A, B), type = "p")
Repita el ejercicio anterior usando como distribución de referencia la distribución uniforme (escriba en el pdf la función creada y la gráfica pedida).
g.betabinom = function(m, n, A, B){ # m: tamano de muestra
x = 0:n
c. = max(betabinom(x, size = n, A, B)/dunif(x, 0, n))
x.m = numeric(0) ; i = 0
while(length(x.m) < m){
j = sample(n+1, size = 1)-1 # paso 1
u = runif(1) # paso 2
if(u <= betabinom(j, size = n, A, B)/(c.*dunif(j, 0, n))){
i = i + 1
x.m[i] = j # paso 3
}
}
x.m
}
# Parametros para la grafica
n = 25
A = 2
B = 5
X = g.betabinom(m = 10000, n, A, B)
plot(table(X)/sum(table(X)))
lines(0:n, betabinom(0:n, n, A, B), type = "p")
El método binomial no parece ser la mejor opción, porque la distribución de masa binomial tiende a cero con valores cercanos a \(x=25\), a psear que estas dos distribuciones tienen forma similar. De esta manera se consigue un valor \(c=46\) cuando \(p=0.432\) (se optimizó \(p\) para obtener el \(c\) más bajo posible). Mientras que con la distribución uniforme se consigue \(c=2.16\), mucho mejor!
Considere la variable aleatoria continua \(X\) con función de densidad de probabilidad
\[ f_{X}(x)=\frac{1}{2 \sigma}\left[1+\cos \left(\frac{x-\mu}{\sigma} \pi\right)\right], \quad \mu-\sigma \leq x \leq \mu+\sigma \] con \(\mu \in \mathbb{R}\) y \(\sigma>0\)
Cree una función que emplee el método de aceptación y rechazo, y use una distribución normal como distribución de referencia, para generar una muestra de tamaño \(m=10000\) de esta distribución (con \(\mu=10\) y \(\sigma=2)\) y grafique además su histograma junto con la función de densidad de probabilidad teórica (escriba en el pdf la función creada y la gráfica pedida).
dcos = function(x, mu, sigma){
1/(2*sigma)*(1+cospi((x-mu)/sigma))
}
rcos = function(m, mu, sigma){ # m: tamano de muestra
x = seq(mu-sigma, mu+sigma, length=m)
c. = max(dcos(x,mu,sigma)/dnorm(x, mu, sigma/3))
x.m = numeric(0) ; i = 0
while(length(x.m) < m){
j = rnorm(1, mu, sigma/3) # paso 1
u = runif(1) # paso 2
if(u <= dcos(j, mu, sigma)/(c.*dnorm(j, mu, sigma/3))){
i = i + 1
x.m[i] = j # paso 3
}
}
x.m
}
mu = 10
sigma = 2
m = 10000
x = seq(mu-sigma, mu+sigma, length=m)
hist(rcos(m, mu, sigma), freq = F, breaks = 50)
curve(dcos(x,mu,sigma), from = min(x), to = max(x), add = T)
Repita el ejercicio anterior usando como distribución de referencia la distribución uniforme (escriba en el pdf la función creada y la gráfica pedida).
rcos = function(m, mu, sigma){ # m: tamano de muestra
x = seq(mu-sigma, mu+sigma, length=m)
c. = max(dcos(x,mu,sigma)/dunif(x, min(x), max(x)))
x.m = numeric(0) ; i = 0
while(length(x.m) < m){
j = runif(1, min(x), max(x)) # paso 1
u = runif(1) # paso 2
if(u <= dcos(j, mu, sigma)/(c.*dunif(j, min(x), max(x)))){
i = i + 1
x.m[i] = j # paso 3
}
}
x.m
}
mu = 10
sigma = 2
m = 10000
x = seq(mu-sigma, mu+sigma, length=m)
hist(rcos(m, mu, sigma), breaks=50, freq = F)
curve(dcos(x,mu,sigma), from = min(x), to = max(x), add = T)
Es más eficiente el método usando como función de referencia la distribución normal, con esta se consigue \(c=1.56\), mientras con la distribución uniforme se consigue \(c=2.33\).
Considere una v.a. con función de densidad de probabilidad
\[f(x)=\frac{\lambda e^{-\lambda x}}{1-e^{-2 \pi \lambda}}, \quad 0 \leq x<2 \pi, \quad \lambda>0\]
Si la función de distribución acumulada de esta v.a. es \[ F(x)=\frac{1-e^{-\lambda x}}{1-e^{-2 \pi \lambda}} \]
construya una función que mediante el método de la transformada inversa genere una muestra de 10000 valores de esta distribución con parámetro \(\lambda=2\) Además, grafique un histograma y la f.d.p. teórica para la secuencia de valores generados. No es necesario agregar los cálculos hechos a mano.
# f.d.p. de X
df = function(x, lambda) {
(lambda*exp(-lambda*x))/(1-exp(-2*pi*lambda))
}
rf = function(u, lambda) {
-log(1-u*(1-exp(-2*pi*lambda)))/lambda
}
n = 10000
lambda = 2
u = runif(n)
x = rf(u, lambda)
hist(x, breaks=50, freq = F)
curve(df(x,lambda), from = 0, to = 2*pi, add = T)