Universidade Federal do Paraná
Curso de Estatística
CE 089 -
Estatística Computacional II - 2014/2
Prof. Dr. Walmes Marques Zeviani
Tabela de conteúdo
##-----------------------------------------------------------------------------
## Amostrador independente.
## Gerar números de uma distribuição Beta usando a distribuição Uniforme
## e/ou normal.
## Distribuição alvo: X ~ Beta.
fX <- function(x) dbeta(x, shape1=2, shape2=3)
curve(fX, 0, 1)
## Distribuição candidata (proposal): X ~ Uniforme.
qX <- function(x) dunif(x, 0, 1)
## Gráfico das densidados sobrepostas.
curve(fX, 0, 1)
curve(qX, add=TRUE, col=2)
legend("topright", legend=c("Alvo", "Candidata"), lty=1, col=1:2, bty="n")
##-----------------------------------------------------------------------------
## Passo a passo.
## Vetor com elementos n elementos.
n <- 10
x <- vector(mode="numeric", length=n)
i <- 2
while(i<n){
## 0. Definir um valor inicial dentro do suporte da distribuição
## alvo.
x[i-1] <- 0.5
## 1. Gerar um número da distribuição candidata.
y <- runif(1); y
## 2. Calcular a quantidade r =
## (fX(y)/fX(x_{i-1}))/(qX(y)/qX(x_{i-1})).
r <- (fX(y)/fX(x[i-1]))/(qX(y)/qX(x[i-1])); r
## 3. Gerar um número aleatório da distribuição uniforme.
u <- runif(1); u
## 4. Se u<r, aceitar o canditado, se não, rejeitar.
if(u<r){
x[i] <- y
print("u<r, então novo valor y na cadeia.")
} else {
x[i] <- x[i-1]
print("u>=r, então último valor da cadeia repetido.")
}
i <- i+1
}
## [1] "u<r, então novo valor y na cadeia."
## [1] "u>=r, então último valor da cadeia repetido."
## [1] "u<r, então novo valor y na cadeia."
## [1] "u<r, então novo valor y na cadeia."
## [1] "u<r, então novo valor y na cadeia."
## [1] "u>=r, então último valor da cadeia repetido."
## [1] "u>=r, então último valor da cadeia repetido."
## [1] "u<r, então novo valor y na cadeia."
##-----------------------------------------------------------------------------
## Fazendo usando gráficos.
## Alvo: Beta.
## Canditada: Uniforme.
iidsampler1 <- function(nsim, x1, plot=FALSE,
go=c("click","enter","none")){
out <- vector(mode="numeric", length=nsim)
## Valor para iniciar a cadeia.
out[1] <- x1
for(i in 2:nsim){
## Realização da distribuição alvo.
if(plot & go[1]=="click"){
y <- locator(n=1)$x
} else {
y <- runif(1)
}
## Cálculo da razão de aceitação.
dg1 <- dbeta(y, 2, 3)
dn1 <- dunif(y)
dg0 <- dbeta(out[i-1], 2, 3)
dn0 <- dunif(out[i-1])
ratio <- (dg1/dg0)/(dn1/dn0)
u <- runif(1)
if(u<ratio){
## Se sim, cadeia ganha novo valor.
out[i] <- y
} else {
## Se não, cadeia recebe o último.
out[i] <- out[i-1]
}
## Parte de representação gráfica do método.
if(plot & nsim<=20){
## Curvas.
curve(dbeta(x, 2, 3), 0, 1, xlim=c(0, 1),
ylab="densidade");
curve(dunif(x), add=TRUE, lty=2);
## Lengendas.
legend("topright",
legend=c(expression(f[X]*" ~ Beta"),
expression(f[Y]*" ~ Unif")),
lty=c(1,2), bty="n")
legend("right",
legend=c(expression("Candidato em"*~i),
expression("Valor em"*~i-1)),
lty=1, col=c(2,4), bty="n")
## Segmentos da base até os valores nas funções.
segments(y, dg1, y, 0, col=2, lty=1);
segments(y, dn1, y, 0, col=2, lty=1);
segments(out[i-1], dg0, out[i-1], 0, col=4, lty=1);
segments(out[i-1], dn0, out[i-1], 0, col=4, lty=1);
## Pontos sobre as funções.
cex <- 2.5; col="yellow"
points(y, dg1, pch=19, cex=cex, col="green");
points(y, dn1, pch=19, cex=cex, col=col);
points(out[i-1], dg0, pch=19, cex=cex, col="green");
points(out[i-1], dn0, pch=19, cex=cex, col=col);
## Rótulos dos pontos.
text(y, dg1, labels=expression(f[X]));
text(y, dn1, labels=expression(f[Y]));
text(out[i-1], dg0, expression(f[X]));
text(out[i-1], dn0, expression(f[Y]));
text(c(y, out[i-1]), 0,
labels=c(expression(x[i]), expression(x[i-1])),
pos=4)
## Anotações matemáticas.
L <- list(dg1=dg1, dg0=dg0, dn1=dn1,
dn0=dn0, num=dg1/dg0, den=dn1/dn0,
ratio=ratio)
L <- lapply(L, round, digits=3)
ex <- substitute(frac(f[X](x[i]), f[X](x[i-1]))/
frac(f[Y](x[i]), f[Y](x[i-1]))*" = "*
frac(dg1, dg0)/frac(dn1, dn0)*" = "*
num/den==ratio, L)
r <- substitute("u = "~u<ratio,
lapply(list(ratio=ratio, u=u),
round, digits=3))
mtext(ex, side=3, line=1, adj=0)
mtext(r, side=3, line=2, adj=1)
mtext(ifelse(u<ratio,
expression(Aceita~x[i]),
expression(Repete~x[i-1])),
side=3, line=1, adj=1)
switch(go[1],
## Avança por cliques do mouse.
click=locator(n=1),
## Avança por enter no console.
console=readline(prompt="Press [enter] to continue"),
## Avança com intervalo de tempo entre etapas.
none=Sys.sleep(0.5))
}
}
return(out)
}
##-----------------------------------------------------------------------------
## Usando.
n <- 10
x <- iidsampler1(n, x1=0.5, plot=TRUE, go="console")
n <- 10
x <- iidsampler1(n, x1=0.5, plot=TRUE, go="click")
##-----------------------------------------------------------------------------
## Alvo: Beta.
## Canditada: Normal.
curve(dbeta(x, 2, 3), 0, 1)
curve(dnorm(x, 0.5, 0.25), add=TRUE, lty=2)
iidsampler2 <- function(nsim, x1, plot=FALSE,
go=c("click","console","none")){
out <- vector(mode="numeric", length=nsim)
out[1] <- x1
for(i in 2:nsim){
## Realização da distribuição alvo.
if(plot & go[1]=="click"){
y <- locator(n=1)$x
} else {
y <- rnorm(1, 0.5, 0.25)
}
## Cálculo da razão de aceitação.
dg1 <- dbeta(y, 2, 3)
dn1 <- dnorm(y, 0.5, 0.25)
dg0 <- dbeta(out[i-1], 2, 3)
dn0 <- dnorm(out[i-1], 0.5, 0.25)
ratio <- (dg1/dg0)/(dn1/dn0)
u <- runif(1)
if(u<ratio){
out[i] <- y
} else {
out[i] <- out[i-1]
}
## Incluir contador da aceitação.
if(plot & nsim<=20){
## Curvas.
curve(dbeta(x, 2, 3), 0, 1, xlim=c(0, 1),
ylab="densidade");
curve(dnorm(x, 0.5, 0.25), add=TRUE, lty=2);
## Lengendas.
legend("topright",
legend=c(expression(f[X]*" ~ Beta"),
expression(f[Y]*" ~ Normal")),
lty=c(1,2), bty="n")
legend("right",
legend=c(expression("Candidato em"*~i),
expression("Valor em"*~i-1)),
lty=1, col=c(2,4), bty="n")
## Segmentos da base até os valores nas funções.
segments(y, dg1, y, 0, col=2, lty=1);
segments(y, dn1, y, 0, col=2, lty=1);
segments(out[i-1], dg0, out[i-1], 0, col=4, lty=1);
segments(out[i-1], dn0, out[i-1], 0, col=4, lty=1);
## Pontos sobre as funções.
cex <- 2.5; col="yellow"
points(y, dg1, pch=19, cex=cex, col="green");
points(y, dn1, pch=19, cex=cex, col=col);
points(out[i-1], dg0, pch=19, cex=cex, col="green");
points(out[i-1], dn0, pch=19, cex=cex, col=col);
## Rótulos dos pontos.
text(y, dg1, labels=expression(f[X]));
text(y, dn1, labels=expression(f[Y]));
text(out[i-1], dg0, expression(f[X]));
text(out[i-1], dn0, expression(f[Y]));
text(c(y, out[i-1]), 0,
labels=c(expression(x[i]), expression(x[i-1])),
pos=4)
## Expressões matemáticas.
L <- list(dg1=dg1, dg0=dg0, dn1=dn1,
dn0=dn0, num=dg1/dg0, den=dn1/dn0,
ratio=ratio)
L <- lapply(L, round, digits=3)
ex <- substitute(frac(f[X](x[i]), f[X](x[i-1]))/
frac(f[Y](x[i]), f[Y](x[i-1]))*" = "*
frac(dg1, dg0)/frac(dn1, dn0)*" = "*
num/den==ratio, L)
r <- substitute("u = "~u<ratio,
lapply(list(ratio=ratio, u=u),
round, digits=3))
mtext(ex, side=3, line=1, adj=0)
mtext(r, side=3, line=2, adj=1)
mtext(ifelse(u<ratio,
expression(Aceita~x[i]),
expression(Repete~x[i-1])),
side=3, line=1, adj=1)
switch(go[1],
click=locator(n=1),
console=readline(prompt="Press [enter] to continue"),
none=Sys.sleep(0.5))
}
}
return(out)
}
##-----------------------------------------------------------------------------
## Usando.
n <- 10
x <- iidsampler2(n, x1=0.5, plot=TRUE, go="none")
##-----------------------------------------------------------------------------
## Gerando muitos números pelo método.
x <- iidsampler2(5000, x1=0.5)
par(mfrow=c(2,2))
plot(x, type="l") ## Traço da cadeia completa.
plot(x[1:100], type="l") ## Traço do começo da cadeia.
acf(x) ## Mostra que a cadeia não é independente.
plot(ecdf(x)) ## Acumulada teórica vs empírica.
curve(pbeta(x, 2, 3), add=TRUE, col=2); layout(1)
##-----------------------------------------------------------------------------
## Alvo: Gama.
## Canditada: Normal.
curve(dgamma(x, 2, 1), -1.5, 8)
curve(dnorm(x, 2/1, sqrt(2/1^2)), add=TRUE, lty=2)
iidsampler3 <- function(nsim, x1, alpha=2, beta=1, mu=NULL, sig=NULL,
plot=FALSE, go=c("click","enter","none")){
out <- vector(mode="numeric", length=nsim)
## Esperança da distribuição proposta.
out[1] <- x1
## Esperança da distribuição proposta.
if(is.null(mu)){
mu <- alpha/beta
}
## Variância da distribuição proposta.
if(is.null(sig)){
sig <- sqrt(alpha/(beta^2))
}
for(i in 2:nsim){
## Realização da distribuição alvo.
if(plot & go[1]=="click"){
y <- locator(n=1)$x
} else {
y <- rnorm(1, mu, sig)
}
## Cálculo da razão de aceitação.
dg1 <- dgamma(y, alpha, beta)
dn1 <- dnorm(y, mu, sig)
dg0 <- dgamma(out[i-1], alpha, beta)
dn0 <- dnorm(out[i-1], mu, sig)
ratio <- (dg1/dg0)/(dn1/dn0)
u <- runif(1)
if(u<ratio){
out[i] <- y
} else {
out[i] <- out[i-1]
}
## Incluir contador da aceitação.
if(plot & nsim<=20){
## Curvas.
curve(dgamma(x, alpha, beta), 0, 8, xlim=c(-2, 8),
ylab="densidade");
curve(dnorm(x, mu, sig), add=TRUE, lty=2);
## Lengendas.
legend("topright",
legend=c(expression(f[X]*" ~ Gama"),
expression(f[Y]*" ~ Normal")),
lty=c(1,2), bty="n")
legend("right",
legend=c(expression("Candidato em"*~i),
expression("Valor em"*~i-1)),
lty=1, col=c(2,4), bty="n")
## Segmentos da base até os valores nas funções.
segments(y, dg1, y, 0, col=2, lty=1);
segments(y, dn1, y, 0, col=2, lty=1);
segments(out[i-1], dg0, out[i-1], 0, col=4, lty=1);
segments(out[i-1], dn0, out[i-1], 0, col=4, lty=1);
## Pontos sobre as funções.
cex <- 2.5; col="yellow"
points(y, dg1, pch=19, cex=cex, col="green");
points(y, dn1, pch=19, cex=cex, col=col);
points(out[i-1], dg0, pch=19, cex=cex, col="green");
points(out[i-1], dn0, pch=19, cex=cex, col=col);
## Rótulos dos pontos.
text(y, dg1, labels=expression(f[X]));
text(y, dn1, labels=expression(f[Y]));
text(out[i-1], dg0, expression(f[X]));
text(out[i-1], dn0, expression(f[Y]));
text(c(y, out[i-1]), 0,
labels=c(expression(x[i]), expression(x[i-1])),
pos=4)
## Expressões matemáticas.
L <- list(dg1=dg1, dg0=dg0, dn1=dn1,
dn0=dn0, num=dg1/dg0, den=dn1/dn0,
ratio=ratio)
L <- lapply(L, round, digits=3)
ex <- substitute(frac(f[X](x[i]), f[X](x[i-1]))/
frac(f[Y](x[i]), f[Y](x[i-1]))*" = "*
frac(dg1, dg0)/frac(dn1, dn0)*" = "*
num/den==ratio, L)
r <- substitute("u = "~u<ratio,
lapply(list(ratio=ratio, u=u),
round, digits=3))
mtext(ex, side=3, line=1, adj=0)
mtext(r, side=3, line=2, adj=1)
mtext(ifelse(u<ratio,
expression(Aceita~x[i]),
expression(Repete~x[i-1])),
side=3, line=1, adj=1)
switch(go[1],
click=locator(n=1),
console=readline(prompt="Press [enter] to continue"),
none=Sys.sleep(0.5))
}
}
return(out)
}
##-----------------------------------------------------------------------------
## Usando.
n <- 10
x <- iidsampler3(n, x1=0.5, plot=TRUE, go="none")
##-----------------------------------------------------------------------------
## Gerando muitos números pelo método.
x <- iidsampler3(5000, x1=0.5)
par(mfrow=c(2,2))
plot(x, type="l") ## Traço da cadeia completa.
plot(x[1:100], type="l") ## Traço do começo da cadeia.
acf(x) ## Mostra que a cadeia não é independente.
plot(ecdf(x)) ## Acumulada teórica vs empírica.
curve(pgamma(x, 2, 1), add=TRUE, col=2); layout(1)
##-----------------------------------------------------------------------------
## Início da cadeia mal escolhido. Não convergência.
set.seed(123)
x <- iidsampler3(5000, x1=9.5)
par(mfrow=c(2,2))
plot(x, type="l") ## Traço da cadeia completa.
plot(x[1:100], type="l") ## Traço do começo da cadeia.
acf(x) ## Mostra que a cadeia não é independente.
plot(ecdf(x)) ## Acumulada teórica vs empírica.
curve(pgamma(x, 2, 1), add=TRUE, col=2); layout(1)
##-----------------------------------------------------------------------------
## Densidade da candidata mal posicionada posicionada.
n <- 5000; alpha <- 2; beta <- 1
mu <- -2; sig <- 2
curve(dgamma(x, alpha, beta), 0, 12, xlim=c(-8, 8), ylab="densidade");
curve(dnorm(x, mu, sig), add=TRUE, lty=2);
set.seed(123)
x <- iidsampler3(n, alpha, beta, mu=mu, sig=sig, x1=2)
par(mfrow=c(2,2))
plot(x, type="l") ## Traço da cadeia completa.
plot(x[1:100], type="l") ## Traço do começo da cadeia.
acf(x) ## Mostra que a cadeia não é independente.
plot(ecdf(x)) ## Acumulada teórica vs empírica.
curve(pgamma(x, alpha, beta), add=TRUE, col=2); layout(1)
##-----------------------------------------------------------------------------
## Metropolis Random Walk.
## Ilustrando o procedimento (com um exemplo bem simples). Obter
## realizações de uma distribuição Normal(0, 1). Definir a distribuição
## candidata qX(x,.) como sendo a uniforme(-delta, delta) que é
## simétrica.
qX <- function(delta, xi){
## delta e xi: escalares parâmetros da distribuição candidata.
## A distribuição candidata é a uniforme.
## Retorna uma realização da distribuição candidata.
runif(1, xi-delta, xi+delta)
}
rwsampler1 <- function(nsim, x1, delta, mu, sigma,
plot=FALSE, go=c("click","enter","none")){
out <- vector(mode="numeric", length=nsim)
out[1] <- x1
for(i in 2:nsim){
## Realização da distribuição alvo.
if(plot & go[1]=="click"){
can <- locator(n=1)$x
} else {
can <- qX(delta, xi=out[i-1])
}
dn1 <- dnorm(can, mu, sigma)
dn0 <- dnorm(out[i-1], mu, sigma)
ratio <- dn1/dn0
u <- runif(1)
if(u<ratio) out[i] <- can else out[i] <- out[i-1]
if(plot & nsim<=20){
curve(dnorm(x, mu, sigma), mu-4*sigma, mu+4*sigma,
ylab="densidade")
curve(dunif(x, out[i-1]-delta, out[i-1]+delta), add=TRUE, lty=2)
du <- dunif(can, out[i-1]-delta, out[i-1]+delta)
## segments(can, du, can, 0, col=4)
segments(can, dn1, can, 0, col=2);
segments(out[i-1], dn0, out[i-1], 0, col=4);
cex <- 2.5; col="yellow"
points(can, dn1, pch=19, cex=cex, col="green");
points(out[i-1], dn0, pch=19, cex=cex, col=col);
## points(can, dn1, pch="N");
## points(out[i-1], dn0, pch="n");
text(can, dn1, expression(f[X]));
text(out[i-1], dn0, expression(f[X]));
ex <- substitute(frac(f[X](x[i]),
f[X](x[i-1]))*" = "*
frac(dn1, dn0)==ratio,
list(dn1=dn1, dn0=dn0, ratio=ratio))
r <- substitute("u = "~u<ratio,
list(ratio=ratio, u=u))
mtext(ex, side=3, line=1, adj=0)
mtext(r, side=3, line=2, adj=1)
mtext(sprintf("então %s", ifelse(u<ratio, "aceita", "rejeita")),
side=3, line=1, adj=1)
switch(go[1],
click=locator(n=1),
console=readline(prompt="Press [enter] to continue"),
none=Sys.sleep(0.5))
}
}
return(out)
}
##-----------------------------------------------------------------------------
## Usando.
mu <- 0; sigma <- 1
x <- rwsampler1(nsim=10, x1=-1, delta=2, mu, sigma, plot=TRUE, go="none")
##-----------------------------------------------------------------------------
## Com muitos valores.
mu <- 0; sigma <- 1
x <- rwsampler1(nsim=5000, x1=0, delta=5, mu, sigma)
par(mfrow=c(2,2))
plot(x, type="l") ## Traço da cadeia completa.
plot(x[1:100], type="l") ## Traço do começo da cadeia.
acf(x) ## Mostra que a cadeia não é independente.
plot(ecdf(x)) ## Acumulada teórica vs empírica.
curve(pnorm(x), add=TRUE, col=2); layout(1)
##-----------------------------------------------------------------------------
## Simular de uma mistura de normais. Normais com variância 1 e mistura 1:1.
k <- 0.5
curve(k*dnorm(x, 0, 1)+(1-k)*dnorm(x, 7, 1), -3, 10)
curve(0.1*dunif(x), add=TRUE, col=2, n=1024)
rwsampler2 <- function(nsim, x1, delta,
plot=FALSE, go=c("click","enter","none")){
out <- vector(mode="numeric", length=nsim)
out[1] <- x1
for(i in 2:nsim){
## Realização da distribuição alvo.
if(plot & go[1]=="click"){
can <- locator(n=1)$x
} else {
can <- qX(delta, xi=out[i-1])
}
dn1 <- k*dnorm(can, 0, 1)+(1-k)*dnorm(can, 7, 1)
dn0 <- k*dnorm(out[i-1], 0, 1)+(1-k)*dnorm(out[i-1], 7, 1)
ratio <- dn1/dn0
u <- runif(1)
if(u<ratio) out[i] <- can else out[i] <- out[i-1]
if(plot & nsim<=20){
curve(k*dnorm(x, 0, 1)+(1-k)*dnorm(x, 7, 1), -3, 10,
ylab="densidade")
curve(0.3*dunif(x, out[i-1]-delta, out[i-1]+delta),
add=TRUE, lty=2)
du <- dunif(can, out[i-1]-delta, out[i-1]+delta)
## segments(can, du, can, 0, col=4)
segments(can, dn1, can, 0, col=2);
segments(out[i-1], dn0, out[i-1], 0, col=4);
cex <- 2.5; col="yellow"
points(can, dn1, pch=19, cex=cex, col="green");
points(out[i-1], dn0, pch=19, cex=cex, col=col);
## points(can, dn1, pch="N");
## points(out[i-1], dn0, pch="n");
text(can, dn1, expression(f[X]));
text(out[i-1], dn0, expression(f[X]));
ex <- substitute(frac(f[X](x[i]),
f[X](x[i-1]))*" = "*
frac(dn1, dn0)==ratio,
list(dn1=dn1, dn0=dn0, ratio=ratio))
r <- substitute("u = "~u<ratio,
list(ratio=ratio, u=u))
mtext(ex, side=3, line=1, adj=0)
mtext(r, side=3, line=2, adj=1)
mtext(sprintf("então %s", ifelse(u<ratio, "aceita", "rejeita")),
side=3, line=1, adj=1)
switch(go[1],
click=locator(n=1),
console=readline(prompt="Press [enter] to continue"),
none=Sys.sleep(0.5))
}
}
return(out)
}
##-----------------------------------------------------------------------------
## Usando.
x <- rwsampler2(nsim=20, x1=1, delta=2, plot=TRUE, go="none")
##-----------------------------------------------------------------------------
## Muitos valores.
## Janela estreita.
set.seed(123)
x <- rwsampler2(nsim=20000, x1=1, delta=0.9, plot=FALSE)
par(mfrow=c(2,2))
plot(x, type="l") ## Traço da cadeia completa.
plot(x[1:100], type="l") ## Traço do começo da cadeia.
acf(x) ## Mostra que a cadeia não é independente.
plot(ecdf(x)) ## Acumulada teórica vs empírica.
curve(k*pnorm(x, 0, 1)+(1-k)*pnorm(x, 7, 1), add=TRUE, col=2); layout(1)
prop.table(table(x<3.5))
##
## FALSE TRUE
## 0.8156 0.1844
## Janela larga.
set.seed(123)
x <- rwsampler2(nsim=20000, x1=1, delta=4, plot=FALSE)
par(mfrow=c(2,2))
plot(x, type="l") ## Traço da cadeia completa.
plot(x[1:100], type="l") ## Traço do começo da cadeia.
acf(x) ## Mostra que a cadeia não é independente.
plot(ecdf(x)) ## Acumulada teórica vs empírica.
curve(k*pnorm(x, 0, 1)+(1-k)*pnorm(x, 7, 1), add=TRUE, col=2); layout(1)
prop.table(table(x<3.5))
##
## FALSE TRUE
## 0.519 0.481