Wie man eine Mischung finiter Komponenten aus verschiedenen parametrischen Familien mit JAGS modelliert?

Question

Stellen Sie sich einen zugrunde liegenden Prozess vor, der eine Zahl aus einer Normalverteilung mit der Wahrscheinlichkeit $ \ alpha $ und aus einer gleichmäßigen Verteilung mit der Wahrscheinlichkeit $ 1 - \ alpha $ zeichnet. Die beobachtete Folge von Zahlen, die durch dieses Verfahren erzeugt folgt daher, eine Verteilung f $ $, die eine Mischung von 2 Komponenten und Mischungsgewichten von $ \ alpha $ und 1 $ - \ alpha $. Wie würden Sie diese Art der Mischung mit JAGS modellieren, wenn die beobachtete Sequenz die einzige Eingabe ist, aber die parametrischen Familien bekannt sind?

Beispiel (in R):

set.seed(8361299) 
N <- 100 
alpha <- 0.3 
mu <- 5 
max <- 50 
# Which component to choose from? 
latent_class <- rbinom(N, 1, alpha) 
Y <- ifelse(latent_class, runif(N, min=mu, max=max), rnorm(N, mean=mu)) 

Die erzeugte (beobachtet) Y wie folgt aussieht:

Mit JAGS, sollte es möglich sein, die Mischungsgewichte zu erhalten, sowie die Parameter, der bekannten Komponenten?

Answer 1

Mischmodelle mit der gleichen parametrischen Verteilung sind in JAGS/BUGS ziemlich einfach, aber Mischmodelle mit variierenden parametrischen Antworten (wie Ihre) sind etwas komplizierter. Eine Methode besteht darin, den "Eindeutigkeitstrick" zu verwenden, wobei wir manuell die Wahrscheinlichkeit der Antwort berechnen (indem wir eine der beiden Verteilungen auswählen, wie durch den latenten Teil des Modells spezifiziert) und diese an die (falsche) Antwort eines Bernoulli-Prozesses anpassen jeder Datenpunkt. Zum Beispiel:

# Your data generation: 
set.seed(8361299) 
N <- 100 
alpha <- 0.3 
mu <- 5 
max <- 50 
# Which component to choose from? 
latent_class <- rbinom(N, 1, alpha) 
Y <- ifelse(latent_class, runif(N, min=mu, max=max), rnorm(N, mean=mu)) 

# The model: 
model <- "model{ 

    for(i in 1:N){ 

     # Log density for the normal part: 
     ld_norm[i] <- logdensity.norm(Y[i], mu, tau) 

     # Log density for the uniform part: 
     ld_unif[i] <- logdensity.unif(Y[i], lower, upper) 

     # Select one of these two densities: 
     density[i] <- exp(ld_norm[i]*norm_chosen[i] + ld_unif[i]*(1-norm_chosen[i])) 

     # Generate a likelihood for the MCMC sampler: 
     Ones[i] ~ dbern(density[i]) 

     # The latent class part as usual: 
     norm_chosen[i] ~ dbern(prob) 

    } 

    # Priors: 
    lower ~ dnorm(0, 10^-6) 
    upper ~ dnorm(0, 10^-6) 
    prob ~ dbeta(1,1) 
    mu ~ dnorm(0, 10^-6) 
    tau ~ dgamma(0.01, 0.01) 

    # Specify monitors, data and initial values using runjags: 
    #monitor# lower, upper, prob, mu, tau 
    #data# N, Y, Ones 
    #inits# lower, upper 
}" 

# Run the model using runjags (or use rjags if you prefer!) 
library('runjags') 

lower <- min(Y)-10 
upper <- max(Y)+10 

Ones <- rep(1,N) 

results <- run.jags(model, sample=20000, thin=1) 
results 

plot(results) 

Dies scheint Ihre Parameter ziemlich gut (Ihr Alpha-1-prob) zu erholen, aber achten Sie auf Autokorrelation (und Konvergenz).

Matt

---

EDIT: Da Sie verallgemeinern, um mehr als zwei Verteilungen gefragt, hier ist äquivalent (aber mehr verallgemeinerbare) Code:

# The model: 
model <- "model{ 
    for(i in 1:N){ 
     # Log density for the normal part: 
     ld_comp[i, 1] <- logdensity.norm(Y[i], mu, tau) 
     # Log density for the uniform part: 
     ld_comp[i, 2] <- logdensity.unif(Y[i], lower, upper) 
     # Select one of these two densities and normalise with a Constant: 
     density[i] <- exp(ld_comp[i, component_chosen[i]] - Constant) 
     # Generate a likelihood for the MCMC sampler: 
     Ones[i] ~ dbern(density[i]) 
     # The latent class part using dcat: 
     component_chosen[i] ~ dcat(probs) 
    } 
    # Priors for 2 parameters using a dirichlet distribution: 
    probs ~ ddirch(c(1,1)) 
    lower ~ dnorm(0, 10^-6) 
    upper ~ dnorm(0, 10^-6) 
    mu ~ dnorm(0, 10^-6) 
    tau ~ dgamma(0.01, 0.01) 
    # Specify monitors, data and initial values using runjags: 
    #monitor# lower, upper, probs, mu, tau 
    #data# N, Y, Ones, Constant 
    #inits# lower, upper, mu, tau 
}" 

library('runjags') 

# Initial values to get the chains started: 
lower <- min(Y)-10 
upper <- max(Y)+10 
mu <- 0 
tau <- 0.01 

Ones <- rep(1,N) 

# The constant needs to be big enough to avoid any densities >1 but also small enough to calculate probabilities for observations of 1: 
Constant <- 10 

results <- run.jags(model, sample=10000, thin=1) 
results 

Dieser Code für so viele Distributionen funktionieren wie Sie benötigen, aber exponentiell schlechtere Autokorrelation mit mehr Verteilungen erwarten.