Bayesianische Schätzung von einer Binomialverteilung

Diese Frage ist eine technische Folge dieser Frage .

Ich habe Probleme, das in Raftery (1988) $N$ vorgestellte Modell zu verstehen und zu replizieren : Inferenz für den Binomial- Parameter: ein hierarchischer Bayes-Ansatz in WinBUGS / OpenBUGS / JAGS. Es geht aber nicht nur um Code, es sollte hier auch zum Thema gehören.

Hintergrund

Sei eine Menge von Erfolgszählungen aus einer Binomialverteilung mit unbekanntem und . Ferner gehe ich davon aus, dass einer Poisson-Verteilung mit dem Parameter folgt (wie in der Veröffentlichung beschrieben). Dann hat jedes eine Poisson-Verteilung mit dem Mittelwert . Ich möchte die Prioren in und angeben . $x=(x_{1},\ldots,x_{n})$ $N$ $\theta$ $N$ $\mu$ $x_{i}$ $\lambda = \mu \theta$ $\lambda$ $\theta$

Unter der Annahme, dass ich keine guten Vorkenntnisse über oder ; habe, möchte ich sowohl als auch ; nicht informative Prioritäten zuweisen . Angenommen, meine Prioren sind und . $N$ $\theta$ $\lambda$ $\theta$ $\lambda\sim \mathrm{Gamma}(0.001, 0.001)$ $\theta\sim \mathrm{Uniform}(0, 1)$

Der Autor verwendet ein falsches Vorzeichen von , WinBUGS akzeptiert jedoch keine falschen Vorzeichen. $p(N,\theta)\propto N^{-1}$

Beispiel

In dem Artikel (Seite 226) werden die folgenden Erfolgszahlen der beobachteten Wasserböcke angegeben: . Ich möchte schätzen , die Größe der Bevölkerung. $53, 57, 66, 67, 72$ $N$

So habe ich versucht, das Beispiel in WinBUGS auszuarbeiten ( aktualisiert nach dem Kommentar von @ Stéphane Laurent):

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)

Das Modell ist Sill nicht konvergieren schön nach 500'000 Proben mit 20'000 Burn-in - Proben. Hier ist die Ausgabe eines JAGS-Laufs:

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

Fragen

Klar, mir fehlt etwas, aber ich kann nicht genau sehen, was. Ich denke, meine Formulierung des Modells ist irgendwo falsch. Meine Fragen sind also:

Warum funktioniert mein Modell und seine Implementierung nicht?
Wie könnte das von Raftery (1988) gegebene Modell richtig formuliert und implementiert werden?

Danke für Ihre Hilfe.

bayesian binomial hierarchical-bayesian winbugs COOLSerdash
quelle

Nach dem Papier sollten Sie hinzufügen mu=lambda/thetaund ersetzen n ~ dpois(lambda)durchn ~ dpois(mu)

Stéphane Laurent

@ StéphaneLaurent Danke für den Vorschlag. Ich habe den Code entsprechend geändert. Leider konvergiert das Modell immer noch nicht.

COOLSerdash

Was passiert, wenn Sie ?

N < 72

$N<72$

Sycorax sagt Reinstate Monica

Wenn

, ist die Wahrscheinlichkeit Null, da Ihr Modell davon ausgeht, dass mindestens 72 Wasserböcke vorhanden sind. Ich frage mich, ob das Probleme für den Sampler verursacht.

N < 72

$N<72$

Sycorax sagt Reinstate Monica

Ich glaube nicht, dass das Problem die Konvergenz ist. Ich denke , das Problem ist , dass der Probennehmer schlecht wegen des hohen Korrelationsgrades bei den verschiedenen Ebenen des Modells

niedrig ist , während

niedrige bezogen auf die Gesamtzahl von Iterationen. Ich würde vorschlagen , nur die posteriore Berechnung direkt, beispielsweise über ein Gitter

\hat{R}

$\hat{R}$

n_{e f f}

$n_{eff}$

θ, N

$\theta, N$

Sycorax sagt Reinstate Monica

Antworten:

Nun, da Sie Ihren Code zum Laufen gebracht haben, scheint diese Antwort etwas zu spät zu sein. Aber ich habe den Code schon geschrieben, also ...

Für das, was es wert ist, ist dies das gleiche * Modell, zu dem es passt rstan. Es wird auf meinem Consumer-Laptop in 11 Sekunden geschätzt, wodurch eine effektivere Stichprobengröße für unsere interessierenden Parameter in weniger Iterationen erreicht wird. $(N, \theta)$

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))

Beachte, dass ich thetaals 2-Simplex gecastet habe . Dies dient nur der numerischen Stabilität. Die Menge des Interesses ist theta[1]; offensichtlich theta[2]ist überflüssige Information.

* Wie Sie sehen, ist die hintere Zusammenfassung praktisch identisch, und die Heraufstufung von zu einer realen Menge scheint keinen wesentlichen Einfluss auf unsere Schlussfolgerungen zu haben. $N$

$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

$N, \theta$ $\tilde{y}$ $\tilde{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

rstan $N$ $\bar{y}=\theta N$

hinten über einem Gitter

Der folgende Code kann bestätigen, dass unsere Ergebnisse aus Stan sinnvoll sind.

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)\times\left\{N|N\in\mathbb{Z}\land N\ge72)\right\}$

Sycorax sagt Reinstate Monica
quelle

+1 und akzeptiert. Ich bin beeindruckt! Ich habe auch versucht, Stan für einen Vergleich zu verwenden, konnte das Modell jedoch nicht übertragen. Die Schätzung meines Modells dauert ungefähr 2 Minuten.

COOLSerdash

Die einzige Schwierigkeit bei Stan für dieses Problem ist, dass alle Parameter real sein müssen, was es ein wenig unpraktisch macht. Aber da Sie die Log-Wahrscheinlichkeit durch eine beliebige Funktion bestrafen können, müssen Sie nur die Mühe machen, sie zu programmieren ... Und die zusammengesetzten Funktionen ausgraben, um dies zu tun ...

Sycorax sagt Reinstate Monica

Ja! Das war genau mein Problem. nkann nicht als Ganzzahl deklariert werden und ich kannte keine Problemumgehung für das Problem.

COOLSerdash

Etwa 2 Minuten auf meinem Desktop.

COOLSerdash

@COOLSerdash Sie könnten an [dieser] [1] Frage interessiert sein, bei der ich frage, welche der Rasterergebnisse oder rstanErgebnisse korrekter sind. [1] stats.stackexchange.com/questions/114366/…

Sycorax sagt Reinstate Monica

$\lambda$

Hier ist mein Analyseskript und die Ergebnisse mit JAGS und 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")

Die Berechnung auf meinem Desktop-PC dauerte ungefähr 98 Sekunden.

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

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

Die Ergebnisse sind:

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

$N$ $522$ $233$ $N$

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

$N$ $150$ $(78; 578)$ $(80; 598)$

COOLSerdash
quelle