O gerador básico de números aleatórios é o que produz valores de uma distribução uniforme, \(U \sum \text{U}[0, 1]\), que vamos cnsiderar disponível.
A partir da geração de valores de \(u\) é possível gerar valores de variáveis com outras distribuições, por exemplo, utilizando o método da transformação inversa. Entretanto, o uso deste método é limitado e outros métodos são necessários e/ou convenientes. No que se segue vamos considerar métodos que exploram relações entre variáveis aleatórias. O material a seguir é baseado em [@rizzo].
Diversar variáveis aleatórias podem ser obtidas por transformações (além da inversa) e relações entre distribuições de variáveis aleatórias. Tais relações sugerem formas de obter simulações de variáveis.
Alguns exemplos:
\[ \begin{align} U,V \stackrel{ind}{\sim} {\rm U}(0,1) \\ & \longrightarrow Z_1 = \sqrt{-2 \log(U)} \; \cos(2 \pi V)\\ & \longrightarrow Z_2 = \sqrt{-2 \log(U)} \; \sin(2 \pi V)\\ Z \sim {\rm N}(0,1) & \longrightarrow V = Z^2 \sim \chi^2(1) \\ Z \sim {\rm N}(0,1) \mbox{ e } V \sim \chi^2(n) & \longrightarrow T = \frac{Z}{\sqrt{V/n}} \sim {\rm t}(n) \\ U \sim \chi^2(m) \mbox{ e } V \sim \chi^2(n) & \longrightarrow F = \frac{U/m}{V/n} \sim {\rm F}(m, n)\\ U \sim {\rm G}(a, \lambda) \mbox{ e } V \sim {\rm G}(b, \lambda) (U, V \; ind.) &\longrightarrow \frac{U}{U+V} \sim {\rm B}(a,b) \end{align} \]
Atividade:
Partindo de um simulações de uniformes (runif()) gerar amostras das distribuições indicadas acima.
Fazer gráficos e obter medidas para avaliar as simulações.
Gerando de \(Y \sim \chi^2(\nu = 1)\). Em R poderíamos uusar diretamente a função () mas, para ilustrar como relações entre variáveis podem ser usadas para definir algorítmos, vamos aqui gerar a partir de valores que seguem uma distribuição normal padrão \(Z \sim \text{N}(0, 1)\). Usamos o resultado:
Valores de estatísticas das amostras devem se aproximar dos teóricos.
## [1] 1.028869
## [1] 2.101801
## [1] 2.101591
Diferentes tipos de gráficos podem ser usados para descrever o comportamento dos dados gerados e comparar com o experado.
par(mfrow=c(1,4), mar=c(3,3,0,0), mgp=c(2,1,0))
hist(y, freq = FALSE, main="")
curve(dchisq(x, df=1), from=0, to=25, col=4, add=T)
##
plot(density(y), main="", xlab = "y")
curve(dchisq(x, df=1), from=0, to=25, col=4, add=T)
#
plot(ecdf(y), main="")
curve(pchisq(x, df=1), from=0, to=25, col=4, add=T)
#
qTEO <- qchisq(ppoints(n), 1)
qqplot(qTEO, y, xlab="quantil teórico")
abline(0, 1)
Certas variáveis podem ser obtidas como soma de outras.
A distribuição da soma de variáveis aleatórias dependentes é a convolução de sua distribuição. A família de distribuições é fechada sob convolução se a distribuição resultante da soma é da mesma família das variáveis originais
Alguns exemplos de distribuições de somas: \[ \begin{align} X \sim {\rm Ber}(p) & \longrightarrow \sum_{i=1}^{n} X_i \sim {\rm Bin}(n, p)\\ X \sim {\rm Geo}(p) & \longrightarrow \sum_{i=1}^{k} X_i \sim {\rm NegBin}(k, p)\\ Y \sim {\rm Exp}(\lambda) & \longrightarrow \sum_{i=1}^{n} X_i \sim {\rm Gama}(n, p) \text{(Erlang para n inteiro)}\\ Z \sim {\rm N}(0,1) & \longrightarrow \sum_{i=1}^{\nu} Z_i^2 \sim \chi^2(\nu) \\ \end{align} \]
Exemplos de distribuições fechadas sob convolução (sob independência) \[ \begin{align} Y_1 \sim {\rm P}(\lambda_1) \text{ e } Y_2 \sim {\rm P}(\lambda_2) & \longrightarrow Y_1 + Y_2 \sim {\rm P}(\lambda_1 + \lambda_2)\\ Y_1 \sim {\rm N}(\mu_1, \sigma^2_1) \text{ e } Y_2 \sim {\rm N}(\mu_1, \sigma^2_1) & \longrightarrow Y_1 + Y_2 \sim {\rm N}(\mu_1 + \mu_2, \sigma^2_1 + \sigma^2_2) \\ Y \sim \chi^2(1) & \longrightarrow \sum_{i=1}^{\nu} Y_i \sim \chi^2(\nu) \\ Y_1 \sim \chi^2(m) \text{ e } Y_2 \sim \chi^2(n) & \longrightarrow Y_1 + Y_2 \sim \chi^2(m+n)\\ Y_i \sim {\rm G}(\alpha_i,\beta) & \longrightarrow \sum_{i=1}^{n} Y_i \sim {\rm G}(\textstyle\sum_{i=1}^{n}\alpha_i,\beta) \\ Y_1 \sim {\rm Cauchy}(\mu_1, \sigma^2_1) \text{ e } Y_2 \sim {\rm Cauchy}(\mu_1, \sigma^2_1) & \longrightarrow Y_1 + Y_2 \sim {\rm Cauchy}(\mu_1 + \mu_2, \sigma^2_1 + \sigma^2_2) \\ \end{align} \]
Gerando agora de \(Y \sim \chi^2(\nu = 3)\) a partir de \(Z \sim \text{N}(0,1)\).
Usamos os resultados:
Valores de estatísticas das amostras devem se aproximar dos teóricos.
## [1] 3.041619
## [1] 6.158102
## [1] 6.157486
Diferentes tipos de gráficos podem ser usados para descrever o comportamento dos dados gerados e comparar com o experado.
par(mfrow=c(1,4), mar=c(3,3,0,0), mgp=c(2,1,0))
hist(y, freq = FALSE, main="")
curve(dchisq(x, df=nu), from=0, to=25, col=4, add=T)
##
plot(density(y), main="", xlab = "y")
curve(dchisq(x, df=nu), from=0, to=25, col=4, add=T)
#
plot(ecdf(y), main="")
curve(pchisq(x, df=nu), from=0, to=25, col=4, add=T)
#
qTEO <- qchisq(ppoints(n), nu)
qqplot(qTEO, y, xlab="quantil teórico")
abline(0, 1)
Neste exemplo vamos contrastar a distribuição gerada por uma convolução (variável \(Y_3\)) e mistura (variável \(Y_4\)) conforme se segue.
Simulando:
\[ \begin{align}
Y_1 &\sim N(\mu = 175, \sigma^2 = 3^3) \\
Y_2 &\sim N(\mu = 165, \sigma^2 = 3^3) \\
Y_3 &= Y_1 + Y_2 \\
Y_{3a} &= \frac{Y_1 + Y_2}{2} \\
Y_4 &= \begin{cases} Y_1 \mbox{ c/ prob } 0.4 \\ Y_2 \mbox{ c/ prob } 0.6 \end{cases}
\end{align}\]
N <- 1000
y1 <- rnorm(N, 175, 3)
y2 <- rnorm(N, 165, 3)
y3 <- y1 + y2
y3a <- (y1 + y2)/2
u <- runif(N)
y4 <- ifelse(u < 0.4, y1, y2)
## forma alternativa de gerar y4
#w <- as.integer(u < 0.4)
#y4 <- w * y1 + (1-w) * y2
As densidades teóricas e empíricas de \(Y_3\), \(Y_{3a}\) e \(Y_4\) são mostradas na figura a seguir. Note que para \(Y_4\) foi necessário definir uma função para calcular a função de distribuição de probabilidades.
par(mfrow=c(1,3), mar=c(3,3,0,0), mgp=c(2,1,0))
plot(density(y3), main="", xlab = "y3")
curve(dnorm(x, m = 340, sd = sqrt(18)), from = 330, to = 360, add = TRUE, col = 4, lwd = 2)
##
plot(density(y3a), main="", xlab = "y3a", xlim = c(155, 185))
curve(dnorm(x, m = 170, sd = sqrt(18/4)), from = 155, to = 185, add = TRUE, col = 4, lwd = 2)
lines(density(y1), lty = 2)
curve(dnorm(x, m = 175, sd = 3), from = 155, to = 185, add = TRUE, col = 4, lty = 2)
lines(density(y2), lty = 2)
curve(dnorm(x, m = 165, sd = 3), from = 155, to = 185, add = TRUE, col = 4, lty = 2)
legend("topright", c("teórica","simulada"), col = c(4, 1), lty = 1)
##
plot(density(y4), main="", xlab = "y3a", xlim = c(155, 185))
mix.f <- Vectorize(function(x, medias, sds, ws)
sum(dnorm(x, m = medias, sd = sds)*ws), "x")
curve(mix.f(x, medias = c(175, 165), sds = c(3,3), ws = c(0.4, 0.6)), add = TRUE, col = 4)
Colocando todas em um mesmo gráfico.
par(mfrow=c(1,1))
plot(density(y3a), col=2, main="", xlim=c(145, 195), xlab = "y")
lines(density(y1))
lines(density(y2))
lines(density(y4), col=6)
legend("topright",
c(expression(Y[1]),expression(Y[2]),expression(Y["3a"]),expression(Y[4])),
lty=1, col=c(1,1,2,6))
Atividade:
Partindo de um simulações de uniformes (runif()) gerar amostras das distribuições indicadas acima.
Fazer gráficos e obter medidas para avaliar as simulações.
Simulando:
\[
\begin{align}
X_1 &\sim G(3, 2) \\
X_2 &\sim G(3, 5) \\
X_3 &= X_1 + X_2\\
X_{3a} &= \frac{X_1 + X_2}{2}\\
X_4 &= \begin{cases} X_1 \mbox{ c/ prob } 0,5 \\ X_2 \mbox{ c/ prob }0,5 \end{cases}
\end{align}
\]
N <- 1000
x1 <- rgamma(N, 3, 2)
x2 <- rgamma(N, 3, 5)
x3 <- x1 + x2
x3a <- (x1 + x2)/2
u <- runif(n)
x4 <- ifelse(u < 0.5, x1, x2)
par(mfcol=c(1,3))
hist(x3, prob=TRUE, main = "")
##
curve(dgamma(x, 3, 5), xlim=c(0,5), col = 4)
curve(dgamma(x, 3, 2), add = TRUE, col = 4)
hist(x3a, prob=TRUE, add = TRUE, col = "transparent")
##
curve(dgamma(x, 3, 5), xlim=c(0,5), col = 4)
curve(dgamma(x, 3, 2), add = TRUE, col = 4)
hist(x4, prob=TRUE, add = TRUE, col = "transparent")
\[ \begin{align} X &= \sum_{j=1}^{5} w_j F(x_j) \\ X_j &\sim {\rm G}(shape=3, rate = 1/j) \; j = 1, \ldots, 5 \\ w_j &= j/15 \end{align} \]
Programação (obs evitando usar for())
n <- 5000
ind <- sample(1:5, size=n, replace=TRUE, prob=(1:5)/15)
rates <- 1/ind
x <- rgamma(n, shape=3, rate=rates)
par(mfrow=c(1,1))
hist(x, main = "", xlim=c(0,35), ylim=c(0,.3), xlab="x", col = "transparent", freq = FALSE, breaks = 25)
for (i in 1:5) curve(dgamma(x, 3, 1/i), from = 0, to = 35, n = 300, col = 4, add = TRUE, lty = 2)
mix.f <- Vectorize(function(x, rates, ws)
sum(dgamma(x, sh = 3, rate = rates)*ws), "x")
curve(mix.f(x, rates = 1/(1:5), ws = (1:5)/15), add = TRUE, col = 4, lw = 2)
legend("topright", c("componentes","mistura","simulação"),
lty = c(2,1,1), lwd = c(1,2,2), col = c(4,4,1))
Atividade
1. Repita último exemplo com pesos iguais para os componentes de mistura.
2. Construa um examplo semelante para mistura de distribuições Beta.
\[ \begin{align} Y &= \sum_{j=1}^{5} w_j F(x_j) \\ Y_j &\sim {\rm G}(shape=3, rate = \beta_j) \; j = 1, \ldots, 5 \\ \beta &= c(1.0, 1.5, 2.0, 2.5, 3.0)^{\prime} w &= c(0.1, 0.2, 0.2, 0.3, 0.2)^{\prime} \end{align} \]
n <- 5000
w <- c(.1,.2,.2,.3,.2)
beta <- c(1,1.5,2,2.5,3)
ind <- sample(1:5, size=n, replace=TRUE, prob=w)
rates <- beta[ind]
x <- rgamma(n, shape=3, rate=rates)
Atividade
Fazer gráfico semelhante ao do exemplo anterior mostrando as componentes, distribuição da mistura e empírica (simulada).
A mistura Poisson-Gama definida a seguir é um examplo de mistura contínua.
\[ \begin{align} Y \sim \text{P}(\lambda) \\ \lambda \sim {\rm G}(3, 2) \end{align} \]
O código a seguir simula desta mistura, ou seja da variável aleatória \(Y\).
alpha <- 3
beta <- 2
N <- 3000
lambda <- rgamma(N, shape = alpha, rate = beta)
y <- rpois(N, lambda)
A distribuição (marginal) \([Y]\) pode ser obtida por: \[ [Y] = \int [Y|\lambda][\lambda] {\rm d}\lambda \] que, neste caso, tem solução analítica levando à distribuição binomial negativa, \[ [Y] \sim \text{NegBin}(\alpha, \frac{\beta}{1+\beta}).\]
Na tabela a seguir comparamos a distribuição de \(Y\) gerada pelo algorítmo baseado na mistura com a distribuição binomial negativa de \(\alpha = 3\) e \(p = \beta/(1+\beta) = 2/3\). As probabilidades obtidas por simulação são proporções estimadas com respectivos erros padrão calcaldos em se.
maxY <- max(y)
mix <- prop.table(table(y))
negbin <- round(dnbinom(0:maxY, alpha, beta/(1+beta)), 3)
se <- sqrt(negbin * (1 - negbin) / N)
round(rbind(mix, negbin, se), 3)
## 0 1 2 3 4 5 6 7 8 9 10
## mix 0.293 0.294 0.197 0.117 0.055 0.026 0.012 0.006 0.001 0.000 0.001
## negbin 0.296 0.296 0.198 0.110 0.055 0.026 0.011 0.005 0.002 0.001 0.000
## se 0.008 0.008 0.007 0.006 0.004 0.003 0.002 0.001 0.001 0.001 0.000
Na figura a segir compara-se a distribuição de \(Y\) gerada pelo algorítmo baseado na mistura com a distribuição Poisson de parâmetro fixo \(\lambda = 1.5\) e com a binomial negativa de \(\alpha = 3\) e \(p = \beta/(1+\beta) = 2/3\).
plot(prop.table(table(y)), type="h", ylim = c(0, 0.35))
lines((0:maxY)+0.1, dpois(0:maxY, lambda=1.5), type="h", col=2)
lines((0:maxY)-0.1, dnbinom(0:maxY, alpha, beta/(1+beta)), type="h", col=4)
legend("topright", c("simulados", "Poisson", "Bin.Neg."), col=c(1,2, 4), lty=1)
Nota-se que a mistura possui valores mais variáveis do que a Poisson com parâmetro fixado e é usual dizer que os valores são superdispersos em relação à distribuição de Poisson.
Misturas deste tipo são usuais na esecificação de modelos de efeitos aleatórios e modelagem Bayesiana.