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

exemplo_gibbs.Rmd

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

```{r}
require(ggplot2)
require(dplyr)
```

```{r}
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
```{r}
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)
}))

```

```{r}
ggplot(as.data.frame(gibbs), aes(x=V1, y=V2)) +
  geom_point() +
  theme_bw()

gibbs %.%
  as.data.frame() %.%
  mutate(iter=1:n()) %.%
  select(phi=V1, m=V2, iter) %.%
  filter(iter > 9000) %.%
  group_by(m) %.%
  mutate(nm=n()) %.%
  ungroup() %.%
  ggplot(aes(x=iter, y=m)) +
  geom_point() +
  geom_step() +
  scale_y_continuous(breaks=1:100 * 2) +
  geom_hline(aes(yintercept=last(m, order_by=nm)), colour='red') +
  theme_bw() +
  theme(axis.text=element_text(size=13))

gibbs %.%
  as.data.frame() %.%
  mutate(iter=1:n()) %.%
  select(phi=V1, m=V2, iter) %.%
  filter(iter > 9000) %.%
  ggplot(aes(x=iter, y=phi)) +
  geom_point() +
  geom_step() +
  scale_y_continuous(breaks=0:20/10) +
  geom_hline(aes(yintercept=mean(phi)), colour='red') +
  theme_bw() +
  theme(axis.text=element_text(size=13))

dados <- gibbs %.%
  as.data.frame() %.%
  mutate(iter=1:n()) %.%
  select(phi=V1, m=V2, iter)
dados[rep(dados$iter, each=2), ] %.%
  mutate(m=lead(m)) %.%
  filter(iter > 9000) %.%
  ggplot(aes(x=phi, y=m)) +
  geom_path() +
  theme_bw()
```

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