diff --git a/animation/amostrador.R b/animation/amostrador.R
new file mode 100644
index 0000000000000000000000000000000000000000..70f890f10704719cca5f54d9961f370ba1ffa35c
--- /dev/null
+++ b/animation/amostrador.R
@@ -0,0 +1,123 @@
+library(animation)
+
+iidsampler <- function(nsim, x1, plot=FALSE,
+ go=c("click", "console", "sleep", "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){
+ par(mfrow=c(1,2))
+ ## Curvas.
+ curve(dbeta(x, 2, 3), 0, 1, xlim=c(0, 1),
+ ylab="Densidade");
+ curve(dunif(x), add=TRUE, lty=2);
+ title(sub=substitute(i==k, env=list(k=i)))
+ ## 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(y[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](y[i]), f[X](x[i-1]))/
+ frac(f[Y](y[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~y[i]~como~x[i]),
+ expression(Repete~x[i-1]~como~x[i])),
+ side=3, line=1, adj=1)
+ plot(out[1:i], xlim=c(0, nsim+1), ylim=c(0,1), type="o")
+ 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.
+ sleep=Sys.sleep(0.5),
+ none=NULL)
+ }
+ }
+ return(out)
+}
+
+iidsampler(n=20, x1=runif(1), plot=TRUE, go="sleep")
+
+##----------------------------------------------------------------------
+
+saveHTML(iidsampler(n=20, x1=runif(1), plot=TRUE, go="none"),
+ img.name="iidsampler",
+ imgdir="iidsampler",
+ interval=0.5,
+ htmlfile="iidsampler.html",
+ ani.width=800, ani.height=400,
+ verbose=FALSE,
+ autoplay=TRUE,
+ autobrowse=FALSE)
+
+saveGIF(iidsampler(n=20, x1=runif(1), plot=TRUE, go="none"),
+ img.name="iidsampler",
+ interval=0.5,
+ movie.name="iidsampler.gif",
+ ani.width=800, ani.height=400,
+ verbose=FALSE, autobrowse=FALSE)
+
+## Apagar diretórios e arquivos.
+file.remove(c("iidsampler.html", "iidsampler.gif"))
+system("rm -r css js iidsampler")
+
+##----------------------------------------------------------------------
diff --git a/animation/randomWalk.R b/animation/randomWalk.R
new file mode 100644
index 0000000000000000000000000000000000000000..7d1e9f25933abd654e2d2aab3f67da05f549dd72
--- /dev/null
+++ b/animation/randomWalk.R
@@ -0,0 +1,285 @@
+##----------------------------------------------------------------------
+## Metropolis Random Walk.
+
+## Ilustrando o procedimento (com um exemplo bem simples). Obter
+## realizações de uma distribuição Normal(0, 1) X~f. Definir a
+## distribuição candidata Y~g como sendo a uniforme(-delta, delta) que é
+## simétrica.
+
+## Avalia a densidade da distribuição alvo.
+f <- function(x, mu=0, sigma=1) dnorm(x, mu, sigma)
+
+## Simula números aleatórios canditatos.
+g <- 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)
+}
+
+curve(f(x, 0, 1), -4, 4)
+curve(dunif(x, 0-2, 0+2), add=TRUE, col=2)
+
+## Vetor com elementos n elementos.
+n <- 10
+x <- vector(mode="numeric", length=n)
+i <- 2
+
+## 1. Definir um valor inicial dentro do suporte da distribuição alvo
+## X.
+x[i-1] <- 0.5
+
+## 2. Gerar um número da distribuição candidata Y.
+y <- g(delta=1, xi=x[i-1]); y
+
+## 3. Calcular a quantidade r = (f(y)/f(x_{i-1}))
+r <- (f(y)/f(x[i-1])); r
+
+## 4. Gerar um número aleatório da distribuição uniforme.
+u <- runif(1); u
+
+## 5. Se u<r, aceitar o canditado, se não, rejeitar e repetir o último
+## número gerado.
+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.")
+}
+
+## 6. Repetir os passos de 2-6.
+i <- i+1
+
+##----------------------------------------------------------------------
+
+rwsampler <- function(nsim, x1,
+ delta, mu, sigma,
+ 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"){
+ can <- locator(n=1)$x
+ } else {
+ can <- g(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){
+ par(mfrow=c(1,2))
+ 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)
+ plot(out[1:i], xlim=c(0, nsim+1), ylim=c(-4,4), type="o",
+ ylab=expression(x[i]), xlab=expression(i))
+ switch(go[1],
+ click=locator(n=1),
+ console=readline(prompt="Press [enter] to continue"),
+ sleep=Sys.sleep(0.5),
+ none=NULL)
+ }
+ }
+ return(out)
+}
+
+##----------------------------------------------------------------------
+## Usando.
+
+mu <- 0; sigma <- 1
+
+## x <- rwsampler(nsim=10, x1=-1, delta=2,
+## mu, sigma, plot=TRUE, go="click")
+
+x <- rwsampler(nsim=20, x1=-1, delta=2,
+ mu, sigma, plot=TRUE, go="none")
+
+library(animation)
+
+saveHTML(rwsampler(nsim=20, x1=-1, delta=2,
+ mu, sigma, plot=TRUE, go="none"),
+ img.name="rw",
+ imgdir="rw",
+ interval=0.5,
+ htmlfile="rw.html",
+ ani.width=800, ani.height=400,
+ verbose=FALSE,
+ autoplay=TRUE,
+ autobrowse=FALSE)
+
+saveGIF(rwsampler(nsim=20, x1=-1, delta=2,
+ mu, sigma, plot=TRUE, go="none"),
+ img.name="rw",
+ interval=0.5,
+ movie.name="rw.gif",
+ ani.width=800, ani.height=400,
+ verbose=FALSE, autobrowse=FALSE)
+
+## Apagar diretórios e arquivos.
+file.remove(c("rw.gif", "rw.html"))
+system("rm -r css js rw")
+
+##----------------------------------------------------------------------
+## Com muitos valores.
+
+mu <- 0; sigma <- 1
+x <- rwsampler(nsim=5000, x1=0, delta=5, mu, sigma, go="none")
+
+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)
+
+rwsamplerMistura <- 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 <- g(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<=50){
+ par(mfrow=c(1,2))
+ 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)
+ plot(out[1:i], xlim=c(0, nsim+1), ylim=c(-3,10), type="o",
+ ylab=expression(x[i]), xlab=expression(i))
+ switch(go[1],
+ click=locator(n=1),
+ console=readline(prompt="Press [enter] to continue"),
+ sleep=Sys.sleep(0.5),
+ none=NULL)
+ }
+ }
+ return(out)
+}
+
+##----------------------------------------------------------------------
+## Usando.
+
+x <- rwsamplerMistura(nsim=20, x1=1, delta=4, plot=TRUE, go="sleep")
+
+saveHTML(rwsamplerMistura(nsim=50, x1=runif(1, 0, 7), delta=4,
+ plot=TRUE, go="none"),
+ img.name="rwm",
+ imgdir="rwm",
+ interval=0.5,
+ htmlfile="rwm.html",
+ ani.width=800, ani.height=400,
+ verbose=FALSE,
+ autoplay=TRUE,
+ autobrowse=FALSE)
+
+saveGIF(rwsamplerMistura(nsim=50, x1=runif(1, 0, 7), delta=4,
+ plot=TRUE, go="none"),
+ img.name="rwm",
+ interval=0.5,
+ movie.name="rwm.gif",
+ ani.width=800, ani.height=400,
+ verbose=FALSE, autobrowse=FALSE)
+
+## Apagar diretórios e arquivos.
+file.remove(c("rwm.gif", "rwm.html"))
+system("rm -r css js rwm")
+
+##----------------------------------------------------------------------
+## Muitos valores.
+
+## Janela estreita.
+set.seed(123)
+x <- rwsamplerMistura(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))
+
+## Janela larga.
+set.seed(123)
+x <- rwsamplerMistura(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))
+
+##----------------------------------------------------------------------