Bayesowskie oszacowanie

16

To pytanie jest techniczną kontynuacją tego pytania .

Mam problem ze zrozumieniem i replikacją modelu przedstawionego w Raftery (1988): Wnioskowanie o dwumianowym parametrze : hierarchiczne podejście BayesaN. w WinBUGS / OpenBUGS / JAGS. Nie chodzi tylko o kod, więc powinien być tutaj na temat.

tło

Niech będzie zbiorem liczb sukcesów z rozkładu dwumianowego o nieznanych i . Ponadto zakładam, że ma rozkład Poissona z parametrem (jak omówiono w artykule). Następnie każdy ma rozkład Poissona ze średnią . Chcę określić priorytety w kategoriach i .N θ N μ x i λ = μ θ λ θx=(x1,,xn)NθNμxiλ=μθλθ

Zakładając, że nie mam żadnej dobrej wcześniejszej wiedzy na temat lub , chcę przypisać nieinformacyjne priory zarówno i . Powiedzmy, że moimi priorytetami są i .θ X θ X ~ G m m ( 0,001 , 0,001 ) θ ~ U n i m O r m ( 0 , 1 )NθλθλGamma(0.001,0.001)θUniform(0,1)

Autor używa niewłaściwego przed ale WinBUGS nie akceptuje niewłaściwych priorytetów.p(N,θ)N1

Przykład

W artykule (strona 226) podano następujące liczby sukcesów zaobserwowanych kozłów wodnych: . Chcę oszacować , liczebność populacji.53,57,66,67,72N

Oto jak próbowałem wypracować przykład w WinBUGS ( zaktualizowany po komentarzu @ Stéphane Laurenta):

model {

# Likelihood
  for (i in 1:N) {
    x[i] ~ dbin(theta, n)
  }

# Priors

n ~ dpois(mu)
lambda ~ dgamma(0.001, 0.001)
theta ~ dunif(0, 1)
mu <- lambda/theta

}

# Data

list(x = c(53, 57, 66, 67, 72), N = 5)

# Initial values

list(n = 100, lambda = 100, theta  = 0.5)
list(n = 1000, lambda = 1000, theta  = 0.8)
list(n = 5000, lambda = 10, theta  = 0.2)

Model ma parapecie nie zbiegają się ładnie po 500.000 próbek 20'000 Burn-in próbek. Oto wynik uruchomienia JAGS:

Inference for Bugs model at "jags_model_binomial.txt", fit using jags,
 5 chains, each with 5e+05 iterations (first 20000 discarded), n.thin = 5
 n.sims = 480000 iterations saved
         mu.vect  sd.vect   2.5%     25%     50%     75%    97.5%  Rhat  n.eff
lambda    63.081    5.222 53.135  59.609  62.938  66.385   73.856 1.001 480000
mu       542.917 1040.975 91.322 147.231 231.805 462.539 3484.324 1.018    300
n        542.906 1040.762 95.000 147.000 231.000 462.000 3484.000 1.018    300
theta      0.292    0.185  0.018   0.136   0.272   0.428    0.668 1.018    300
deviance  34.907    1.554 33.633  33.859  34.354  35.376   39.213 1.001  43000

pytania

Najwyraźniej czegoś mi brakuje, ale nie widzę, co dokładnie. Myślę, że moje sformułowanie modelu jest gdzieś błędne. Więc moje pytania to:

  • Dlaczego mój model i jego wdrożenie nie działają?
  • Jak poprawnie sformułować i wdrożyć model podany przez Raftery (1988)?

Dzięki za pomoc.

COOLSerdash
źródło
2
Po papierze należy dodać mu=lambda/thetai wymienić n ~ dpois(lambda)zn ~ dpois(mu)
Stéphane Laurent
@ StéphaneLaurent Dzięki za sugestię. Odpowiednio zmieniłem kod. Niestety model nadal nie jest zbieżny.
COOLSerdash
1
Co dzieje się, gdy próbkujesz ? N.<72
Sycorax mówi Przywróć Monikę
1
Jeśli , prawdopodobieństwo wynosi zero, ponieważ twój model zakłada, że ​​jest co najmniej 72 kozioł wodny. Zastanawiam się, czy to powoduje problemy z samplerem. N.<72
Sycorax mówi Przywróć Monikę
3
Nie sądzę, że problemem jest konwergencja. Że problemem jest to, że próbnik jest nieskuteczne z powodu wysokiego stopnia korelacji z wielu poziomów jest niska, a n e f f jest niska w stosunku do całkowitej liczby iteracji. Proponuję po prostu obliczania posterior bezpośrednio, na przykład, przez siatki θ , N . R^neffθ,N
Sycorax mówi Przywróć Monikę

Odpowiedzi:

7

Cóż, skoro masz kod do działania, wygląda na to, że ta odpowiedź jest trochę za późna. Ale już napisałem kod, więc ...

Jeśli chodzi o wartość, jest to ten sam * model, który pasuje rstan. Szacuje się to na 11 sekund na moim laptopie konsumenckim, osiągając wyższy efektywny rozmiar próbki dla naszych interesujących parametrów w mniejszej liczbie iteracji.(N,θ)

raftery.model   <- "
    data{
        int     I;
        int     y[I];
    }
    parameters{
        real<lower=max(y)>  N;
        simplex[2]      theta;
    }
    transformed parameters{
    }
    model{
        vector[I]   Pr_y;

        for(i in 1:I){
            Pr_y[i] <-  binomial_coefficient_log(N, y[i])
                        +multiply_log(y[i],         theta[1])
                        +multiply_log((N-y[i]),     theta[2]);
        }
        increment_log_prob(sum(Pr_y));
        increment_log_prob(-log(N));            
    }
"
raft.data           <- list(y=c(53,57,66,67,72), I=5)
system.time(fit.test    <- stan(model_code=raftery.model, data=raft.data,iter=10))
system.time(fit     <- stan(fit=fit.test, data=raft.data,iter=10000,chains=5))

Zauważ, że rzucam thetajako 2-simpleks. To tylko dla stabilności numerycznej. Kwota odsetek wynosi theta[1]; oczywiście theta[2]jest to zbędna informacja.

* Jak widać, tylne podsumowanie jest praktycznie identyczne, a promowanie do rzeczywistej ilości nie wydaje się mieć istotnego wpływu na nasze wnioski.N

Kwantyl 97,5% dla jest o 50% większy dla mojego modelu, ale myślę, że dzieje się tak, ponieważ próbnik Stana jest lepszy w badaniu pełnego zakresu tylnej niż zwykły chód losowy, więc może łatwiej dostać się w reszkę. Jednak mogę się mylić.N

            mean se_mean       sd   2.5%    25%    50%    75%   97.5% n_eff Rhat
N        1078.75  256.72 15159.79  94.44 148.28 230.61 461.63 4575.49  3487    1
theta[1]    0.29    0.00     0.19   0.01   0.14   0.27   0.42    0.67  2519    1
theta[2]    0.71    0.00     0.19   0.33   0.58   0.73   0.86    0.99  2519    1
lp__      -19.88    0.02     1.11 -22.89 -20.31 -19.54 -19.09  -18.82  3339    1

Biorąc wartości wygenerowane ze stanu, używam ich do narysowania późniejszych wartości predykcyjnych ˜ y . Nie powinniśmy się dziwić, że średnia prognoz późniejszych ˜ y jest bardzo zbliżona do średniej z danych przykładowych!N,θy~y~

N.samples   <- round(extract(fit, "N")[[1]])
theta.samples   <- extract(fit, "theta")[[1]]
y_pred  <- rbinom(50000, size=N.samples, prob=theta.samples[,1])
mean(y_pred)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  32.00   58.00   63.00   63.04   68.00  102.00 

rstanNy¯=θN

z tyłu nad siatką

Poniższy kod może potwierdzić, że nasze wyniki ze stanu mają sens.

theta   <- seq(0+1e-10,1-1e-10, len=1e2)
N       <- round(seq(72, 5e5, len=1e5)); N[2]-N[1]
grid    <- expand.grid(N,theta)
y   <- c(53,57,66,67,72)
raftery.prob    <- function(x, z=y){
    N       <- x[1]
    theta   <- x[2]
    exp(sum(dbinom(z, size=N, prob=theta, log=T)))/N
}

post    <- matrix(apply(grid, 1, raftery.prob), nrow=length(N), ncol=length(theta),byrow=F)    
approx(y=N, x=cumsum(rowSums(post))/sum(rowSums(post)), xout=0.975)
$x
[1] 0.975

$y
[1] 3236.665

rstan(0,1)×{N|NZN72)}

Sycorax mówi Przywróć Monikę
źródło
+1 i zaakceptowane. Jestem pod wrażeniem! Próbowałem też użyć Stana do porównania, ale nie mogłem przenieść modelu. Szacowanie mojego modelu zajmuje około 2 minut.
COOLSerdash
Jedyny problem związany ze stanem tego problemu polega na tym, że wszystkie parametry muszą być rzeczywiste, więc jest to trochę niewygodne. Ale ponieważ możesz ukarać prawdopodobieństwo dziennika za pomocą dowolnej arbitralnej funkcji, po prostu musisz przejść przez kłopot, aby go zaprogramować ... I wykopać złożone funkcje, aby to zrobić ...
Sycorax mówi Przywróć Monikę
Tak! To był dokładnie mój problem. nnie można zadeklarować jako liczbę całkowitą i nie znałem obejścia problemu.
COOLSerdash
Około 2 minut na moim pulpicie.
COOLSerdash
1
@COOLSerdash Być może zainteresuje Cię [to] [1] pytanie, w którym pytam, które wyniki siatki lub rstanwyniki są bardziej poprawne. [1] stats.stackexchange.com/questions/114366/…
Sycorax mówi Przywróć Monikę
3

λ

Oto mój skrypt analizy i wyniki przy użyciu JAGS i R:

#===============================================================================================================
# Load packages
#===============================================================================================================

sapply(c("ggplot2"
         , "rjags"
         , "R2jags"
         , "hdrcde"
         , "runjags"
         , "mcmcplots"
         , "KernSmooth"), library, character.only = TRUE)

#===============================================================================================================
# Model file
#===============================================================================================================

cat("
    model {

    # Likelihood    
    for (i in 1:N) {
      x[i] ~ dbin(theta, n)
    }

    # Prior       
    n ~ dpois(mu)
    lambda ~ dgamma(0.005, 0.005)
#     lambda ~ dunif(0, 1000)
    mu <- lambda/theta
    theta ~ dunif(0, 1)    
}    
", file="jags_model_binomial.txt")


#===============================================================================================================
# Data
#===============================================================================================================

data.list <- list(x = c(53, 57, 66, 67, 72, NA), N = 6) # Waterbuck example from Raftery (1988)

#===============================================================================================================
# Inits
#===============================================================================================================

jags.inits <- function() { 
  list(
    n = sample(max(data.list$x, na.rm = TRUE):1000, size = 1) 
    , theta = runif(1, 0, 1)
    , lambda = runif(1, 1, 10)
#     , cauchy  = runif(1, 1, 1000)
    #     , mu = runif(1, 0, 5)
  )
}

#===============================================================================================================
# Run the chains
#===============================================================================================================

# Parameters to store

params <- c("n"
            , "theta"
            , "lambda"
            , "mu"
            , paste("x[", which(is.na(data.list[["x"]])), "]", sep = "")
)

# MCMC settings

niter <- 500000 # number of iterations
nburn <- 20000  # number of iterations to discard (the burn-in-period)
nchains <- 5    # number of chains

# Run JAGS

out <- jags(
  data                 = data.list
  , parameters.to.save = params
  , model.file         = "jags_model_binomial.txt"
  , n.chains           = nchains
  , n.iter             = niter
  , n.burnin           = nburn
  , n.thin             = 50
  , inits              = jags.inits
  , progress.bar       = "text")

Obliczenia trwały około 98 sekund na moim komputerze stacjonarnym.

#===============================================================================================================
# Inspect results
#===============================================================================================================

print(out
      , digits = 2
      , intervals = c(0.025, 0.1, 0.25, 0.5, 0.75, 0.9,  0.975))

Wyniki są następujące:

Inference for Bugs model at "jags_model_binomial.txt", fit using jags,
 5 chains, each with 5e+05 iterations (first 20000 discarded), n.thin = 50
 n.sims = 48000 iterations saved
         mu.vect sd.vect  2.5%    10%    25%    50%    75%     90%   97.5% Rhat n.eff
lambda     62.90    5.18 53.09  56.47  59.45  62.74  66.19   69.49   73.49    1 48000
mu        521.28  968.41 92.31 113.02 148.00 232.87 467.10 1058.17 3014.82    1  1600
n         521.73  968.54 95.00 114.00 148.00 233.00 467.00 1060.10 3028.00    1  1600
theta       0.29    0.18  0.02   0.06   0.13   0.27   0.42    0.55    0.66    1  1600
x[6]       63.03    7.33 49.00  54.00  58.00  63.00  68.00   72.00   78.00    1 36000
deviance   34.88    1.53 33.63  33.70  33.85  34.34  35.34   36.81   39.07    1 48000

N522233N

jagsfit.mcmc <- as.mcmc(out)
jagsfit.mcmc <- combine.mcmc(jagsfit.mcmc)

hpd.80 <- hdr.den(log(as.vector(jagsfit.mcmc[, "n"])), prob = c(80), den = bkde(log(as.vector(jagsfit.mcmc[, "n"])), gridsize = 10000))

exp(hpd.80$mode)

[1] 149.8161

N

(hpd.ints <- HPDinterval(jagsfit.mcmc, prob = c(0.8)))

               lower      upper
deviance 33.61011007  35.677810
lambda   56.08842502  69.089507
mu       72.42307587 580.027182
n        78.00000000 578.000000
theta     0.01026193   0.465714
x[6]     53.00000000  71.000000

N150(78;578)(80;598)

COOLSerdash
źródło