Exemplo amostrador de gibbs com poisson com ponto de mudança.

exemplo_gibbs.R

# Exemplo de amostrador de gibbs

# Ronaldo, quando era magro, marcava um número de gols por mês que seguia distribuição Poisson(2).
# Ao passar do tempo, teve um momento em que o número de gols de Ronaldo mudou de distribuição.
# Essa distribuição continuou sendo poisson, mas mudou o parâmetro...
# Eu acho que essa mudança pode ter acontecido em qualquer mês entre os anos de 1998 e 2010, de forma
# igual, e também ouso dizer que a média dele agora não deve ser muito diferente de uma exponencial(1).
# pergunto: quando foi que o infeliz (mas rico) mudou a forma de jogar, e pra quanto foi a sua média de gols?
# minha função de perda para acertar o mês é 0-1 e minha função de perda para acertar a média é quadrática.

y <- c(1L, 1L, 2L, 4L, 1L, 4L, 4L, 2L, 2L, 0L, 1L, 1L, 3L, 1L, 3L, 
       2L, 3L, 6L, 1L, 3L, 4L, 1L, 2L, 0L, 1L, 1L, 0L, 1L, 4L, 1L, 2L, 
       2L, 2L, 1L, 3L, 2L, 3L, 0L, 3L, 2L, 3L, 2L, 3L, 2L, 2L, 3L, 0L, 
       2L, 1L, 1L, 1L, 2L, 1L, 0L, 0L, 0L, 0L, 1L, 1L, 1L, 2L, 0L, 1L, 
       0L, 1L, 0L, 1L, 2L, 0L, 2L, 0L, 3L, 1L, 0L, 1L, 2L, 2L, 1L, 2L, 
       3L, 1L, 1L, 1L, 0L, 2L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 2L, 2L, 
       2L, 1L, 1L, 2L, 1L, 1L, 0L, 0L, 4L, 2L, 0L, 0L, 1L, 3L, 1L, 3L, 
       1L, 0L, 1L, 0L, 0L, 1L, 0L, 1L, 1L)
a <- 1
lamb <- 2
n <- 120

# Sabemos as condicionais completas para m e phi
rcond_m <- function(y, lamb, phi, a, n) {
  pm <- sapply(1:n, function(x) lamb^(sum(y[1:x])) * exp(-x*lamb) * phi^(sum(y)-sum(y[1:x]+a-1)) * exp(-phi*(n-x+1)))
  r <- sample(1:n, 1, prob=pm/sum(pm))
}
rcond_phi <- function(y, lamb, m, a, n) {
  rgamma(1, shape=a+sum(y)-sum(y[1:m]), rate=1+n-m)
}

m <- 60 # valor inicial
phi <- NULL
gibbs <- t(replicate(10000,  {
  phi <<- rcond_phi(y, lamb, m, a, n)
  m <<- rcond_m(y, lamb, phi, a, n)
  c(phi, m)
}))

# posteriori (com burn-in de 9000)
mean(gibbs[-c(1:9000),1])
as.numeric(names(sort(table(gibbs[-c(1:9000),2]), decreasing=T)[1]))