#========================================================================================== # ce718 - Métodos Computacionalmente Intensivos Paulo Justiniano # Gibbs Sampler (08/08/2011) # Walmes Zeviani #========================================================================================== curve(dchisq(x, df=3), 0, 10) hi <- function(x, sigma2){ (x/sigma2)*exp(-x^2/(2*sigma2)) } curve(hi(x, sigma2=1), 0, 10) curve(dchisq(x, df=1), add=TRUE, col=2) #------------------------------------------------------------------------------------------ # Gibbs Sampler - Amostrador de Gibbs # *permite gerar realizações de uma distribuição de probabilidade alvo; # *a cadeira de valores gerada converge para a distribuição alvo; # *os valores são gerados das suas distribuições condicionais completas (crítico); #------------------------------------------------------------------------------------------ # ilustração do procedimento (com um exemplo muito ridículo): # obtendo realizações da distribuição normal bivariada N(mu=c(0,0), sd=c(1,1), cor=0.5) # função que gra valores das condicionais completas de X1|X2=x2 e X2|X1=x1 (são iguais) ccX <- function(x, cor){ rnorm(1, cor*x, sqrt(1-cor^2)) } nsim <- 2000 # número de simulações cor <- 0.95 # trocar valor muda correlação da série x <- matrix(0, nrow=nsim, ncol=2) for(i in 2:nsim){ x[i-1,2] <- ccX(x[i-1,1], cor) # atualiza x2(i) usando x1(i) x[i,1] <- ccX(x[i-1,2], cor) # atualiza x1(i+1) usando x2(i) } # fazer essa figura animada # fazer um metropolis bivariado para o mesmo problema plot(x, type="s") points(x, pch=19) text(x, labels=1:nsim, pos=sample(1:4, nsim, repl=TRUE)) plot(x[,1], type="l") # gráfico da série plot(acf(x[,1])) # gráfico da correlação da série apply(x, 2, mean) var(x) # P(X1<0 e X2<0) p1 <- x[,1]<0 p2 <- x[,2]<0 p12 <- p1*p2 plot(cumsum(p12)/1:nsim, type="l") #------------------------------------------------------------------------------------------ # estimação de média e variância para uma amostra com especificação de priores # x é uma amostra aleatória ~N(mu, sigma2) # priores: mu~N(m,s2) e sigma-2=phi~gamma(a,b) # nessa especificação a posteriori conjunta p(mu,sigma2|x) não tem forma fechada # usando o amostrador de Gibbs: # 1) obter uma realização para mu da distribuição p(mu|sigma2,x), é a condicional completa; # 2) obter uma realização para sigma2 da distribuição p(sigma2|mu,x); # as condicionais completas para mu e sigma2 são: # mu: normal(m*, s2*), com m*=(m/s2+n*xbar/sigma2)/(1/s2+1/sigma2) e s2*=1/(1/s2+n/sigma2) # sigma-2=phi: gamma(a*, b*) com a*=a+n/2 e b*=sum(1:n) (xi-mu)^2/2 + b # funções para as condicionais completas ccmu <- function(m, s2, sigma2, xbar, n){ ## depende do valor desconhecido sigma2 ## depende dos valores conhecidos m, s2, xbar, n m_star <- (m/s2+n*xbar/sigma2)/(1/s2+n/sigma2) s2_star <- 1/(1/s2+n/sigma2) ## retorna uma realização de mu de sua condicional completa rnorm(1, m_star, s2_star) } ccsigma2 <- function(a, n, x, mu, b){ ## depende do valor desconhecido mu ## depende dos valores conhecidos a, n, x, b a_star <- a+n/2; b_star <- sum((x-mu)^2)/2+b phi <- rgamma(1, a_star, b_star) ## retorna uma realização de sigma2 de sua condicional completa 1/phi } nsim <- 20000 # número de simulações th <- matrix(c(1,8), nrow=nsim, ncol=2, dimnames=list(NULL, c("sigma2","mu"))) set.seed(123) x <- rnorm(10, 8, sqrt(1.5)) # amostra aleatória n <- length(x); xbar <- mean(x) # estatísticas m <- 10; s2 <- 30; # hiperparametros da priori de mu a <- 7; b <- 4; # hiperparametros da priori de sigma-2=phi curve(dgamma(x,7,4), 0, 5) for(i in 2:nsim){ th[i-1,2] <- ccmu(m=m, s2=s2, sigma2=th[i-1,1], xbar=xbar, n=n) th[i,1] <- ccsigma2(a=a, n=n, x=x, mu=th[i-1,2], b=b) } show <- 1:20 plot(th[show,], type="s") points(th[show,], pch=19) text(th[show,], labels=show, pos=sample(1:4, length(show), repl=TRUE)) plot(th) apply(th, 2, mean) apply(th, 2, sd) apply(th, 2, quantile, p=c(0.025, 0.975)) plot(th[,1], type="l") plot(th[,2], type="l") acf(th[,2]) acf(th[,1]) #------------------------------------------------------------------------------------------ # y: dados de número de falhas para cada uma de 10 bombas em uma usina nuclear # t: tempo em observação de cada uma das bombas y <- c(5, 1, 5, 14, 3, 19, 1, 1, 4, 22) t <- c(94, 16, 63, 126, 5, 31, 1, 1, 2, 10) rbind(y, t) # modelar o número de eventos por unidade de tempo e por bomba # verossimilhança: poisson(lambda_i*t_i) # prioris: lambda_i ~Gamma(alpha, beta), alpha=1.8 # beta ~Gamma(gamma, delta), gamma=0.01, delta=1 # fazer inferências para os 11 parâmetros (10 lambdas, 1 beta) # posteriori: p(lambda, beta|y,t) =~ (prod poisson(lambda_i*t_i)*Gamma(alpha*beta))*Gamma(gamma, delta) lambda <- rep(1, 10) post <- function(nsim, beta.ini, alpha, gamma, delta, y, t){ ## nsim: escalar número de simulações ## beta.ini: escalar valor inicial para beta ## alpha: escalar valor do hiperparâmetro da priori de lambda_i ## gamma: escalar valor do hiperparâmetro da priori de beta ## delta: escalar valor do hiperparâmetro da priori de beta ## y: vetor de valores da resposta observada ## t: tempode observação associado à cada valor beta.serie <- c() # para armazenar a série de beta lambda.serie <- matrix(NA, nrow=nsim, ncol=length(y)) # para armazenar a série de lambda beta.at <- beta.ini cc.lambda <- function(alpha, beta, y, t){ # função que gera uma realização rgamma(length(y), y+alpha, t+beta) # de lambda_i da sua condicional completa } cc.beta <- function(alpha, gamma, delta, lambda, y){ # função que gera uma realização rgamma(1, length(y)*alpha+gamma, delta+sum(lambda)) # de beta da sua condicional completa } for(i in 1:nsim){ lambda.at <- cc.lambda(alpha=alpha, beta=beta.at, y=y, t=t) beta.at <- cc.beta(alpha=alpha, gamma=gamma, delta=delta, lambda=lambda.at, y=y) beta.serie[i] <- beta.at lambda.serie[i,] <- lambda.at } return(list(lambda=lambda.serie, beta=beta.serie)) } post.amostra <- post(10000, beta.ini=1, alpha=1.8, gamma=0.01, delta=1, y=y, t=t) apply(post.amostra$lambda, 2, mean) theta <- data.frame(cbind(post.amostra$lambda, post.amostra$beta)) names(theta) <- c(paste("lambda_", 1:10, sep=""), "beta") require(lattice) idx <- sample(length(post.amostra$beta), 200) splom(theta[idx,]) apply(theta, 2, function(x){ c(m=mean(x), s=sd(x), quantile(x, p=c(25,500,975)/1000)) }) theta <- stack(theta) densityplot(~values|ind, data=theta, scales="free", pch=NA, col=1, panel=function(x, subscripts, ...){ panel.densityplot(x, subscripts=subscripts, ...) panel.abline(v=c(mean(x), quantile(x, p=c(2.5,50,97.5)/100)), lty=c(1,2,2,2)) panel.rug(x, col=1) }) show <- 1:20 plot(post.amostra$lambda[show,1]~post.amostra$beta[show], type="s") points(post.amostra$lambda[show,1]~post.amostra$beta[show], pch=19) text(post.amostra$beta[show], post.amostra$lambda[show,1], label=show, pos=sample(1:4, length(show), repl=TRUE)) acf(post.amostra$lambda[,1]) acf(post.amostra$beta) #------------------------------------------------------------------------------------------ # TO DO: Gibbs sampling for Poisson change point problem in R browseURL("http://www-m4.ma.tum.de/courses/SS09/cs/gibbs-pois.pdf") browseURL("http://www.stat.missouri.edu/~ferreiram/stat9250/changepoint.r") #------------------------------------------------------------------------------------------