Universidade Federal do Paraná
Curso de Estatística
CE 089 - Estatística Computacional II - 2014/2
Prof. Dr. Walmes Marques Zeviani


Aula 06

Tabela de conteúdo



Amostrador independente

##-----------------------------------------------------------------------------
## 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")

plot of chunk unnamed-chunk-2

##-----------------------------------------------------------------------------
## 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)

plot of chunk unnamed-chunk-5

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)

plot of chunk unnamed-chunk-7

##-----------------------------------------------------------------------------
## 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)

plot of chunk unnamed-chunk-8

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")