如何使用 JAGS 对来自不同参数族的有限组件的混合进行建模?
How to model a mixture of finite components from different parametric families with JAGS?
想象一个基础过程,它从概率为 $\alpha$ 的正态分布和概率为 $1 - \alpha$ 的均匀分布中抽取数字。
因此,观察到的由该过程生成的数字序列遵循分布 $f$,它是 2 components 和 mixing 的 mixture $\alpha$ 和 $1 - \alpha$ 的权重。
当观察到的序列是唯一输入但参数族已知时,您将如何使用 JAGS 对这种混合进行建模?
示例(在 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))
生成的(观察到的)Y 看起来像:
有了JAGS,应该可以得到混合权重,以及已知成分的参数?
相同参数分布的混合模型在 JAGS/BUGS 中非常简单,但具有不同参数响应的混合模型(如您的模型)有点棘手。一种方法是使用 'ones trick' 我们手动计算响应的可能性(选择模型的潜在部分指定的两个分布之一)并将其与伯努利试验的(假)响应相匹配对于每个数据点。例如:
# 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)
这似乎很好地恢复了你的参数(你的 alpha 是 1-prob),但要注意自相关(和收敛)。
马特
编辑:既然你询问了关于推广到 2 个以上的分布,这里是等效的(但更可推广的)代码:
# 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
此代码适用于您需要的任意数量的分布,但预计随着更多分布的自相关性呈指数级恶化。
想象一个基础过程,它从概率为 $\alpha$ 的正态分布和概率为 $1 - \alpha$ 的均匀分布中抽取数字。 因此,观察到的由该过程生成的数字序列遵循分布 $f$,它是 2 components 和 mixing 的 mixture $\alpha$ 和 $1 - \alpha$ 的权重。 当观察到的序列是唯一输入但参数族已知时,您将如何使用 JAGS 对这种混合进行建模?
示例(在 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))
生成的(观察到的)Y 看起来像:
有了JAGS,应该可以得到混合权重,以及已知成分的参数?
相同参数分布的混合模型在 JAGS/BUGS 中非常简单,但具有不同参数响应的混合模型(如您的模型)有点棘手。一种方法是使用 'ones trick' 我们手动计算响应的可能性(选择模型的潜在部分指定的两个分布之一)并将其与伯努利试验的(假)响应相匹配对于每个数据点。例如:
# 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)
这似乎很好地恢复了你的参数(你的 alpha 是 1-prob),但要注意自相关(和收敛)。
马特
编辑:既然你询问了关于推广到 2 个以上的分布,这里是等效的(但更可推广的)代码:
# 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
此代码适用于您需要的任意数量的分布,但预计随着更多分布的自相关性呈指数级恶化。