Установка смеси распределений (гауссовы + униформа) в Winbugs

Question

Я пытаюсь подгонять модель распределения смеси к вектору значений, смесь должна состоять из 2 гауссовских распределений и 1 равномерного распределения. Я пытаюсь реализовать это в Winbugs. Я нашел много примеров, которые использовали смесь гауссианцев, но не могут понять, как добавить форму. В приведенном ниже паспорте кода параметр параметризуется, чтобы соответствовать векторам значений, масштабированных между нулем и одним, но я получаю «несколько определений узла NSD [1]», поэтому кажется, что моя структура по-прежнему неправа. Какие-либо предложения?

model{ 

    ## priors 
    xmin~dunif(0,1) 
    eps2~dunif(0,1) 
    xmax<-min(xmin+eps2, 1) 
    mu1~dunif(0,1) 
    eps1~dunif(0,1) 
    mu2<-min(mu1+eps1,1) 

    sigma1 ~ dunif(0,.5)  
    sigma2 ~ dunif(0,.5)  
    tau1<-pow(sigma1,-2) 
    tau2<-pow(sigma2,-2) 
    alpha[1]<-1 
    alpha[2]<-1 
    alpha[3]<-1 
    p.state[1:3]~ddirch(alpha[]) 

    for (t in 1:npts) { 
    idx[t] ~ dcat(p.state[]) ## idx is the latent variable and the parameter index 
    x[t,1]~dnorm(mu1,tau1) 
    x[t,2]~dnorm(mu2,tau2) 
    x[t,3]~dunif(xmin,xmax) 

     NSD[t] <-x[t,idx[t]]  
     } 
} 

Answer 1

Вы можете попробовать использовать неинформативные dnorm до вместо dunif предварительного, так что вы можете моделировать до для NSD как ~ dnorm(mu[idx[t]], tau[idx[t]]). Тем не менее, вам нужно усечь, так что можно установить очень низкие/высокие границы для усечения, когда вы хотите нормальных приоритетов.

Может быть что-то вроде этого:

model { 
    mu[1] ~ dunif(0, 1) 
    mu[2] <- min(mu[1] + eps[1], 1) 
    mu[3] <- 0.5 
    eps[1] ~ dunif(0, 1) 
    eps[2] ~ dunif(0, 1) 
    sigma[1] ~ dunif(0,.5)  
    sigma[2] ~ dunif(0,.5)  
    tau[1] <- pow(sigma[1],-2) 
    tau[2] <- pow(sigma[2],-2) 
    tau[3] <- 0.000001 
    left[1] <- -100 # something relatively very low 
    left[2] <- -100 # something relatively very low 
    left[3] ~ dunif(0, 1) 
    right[1] <- 100 # something relatively very high 
    right[2] <- 100 # something relatively very high 
    right[3] <- min(left[3] + eps[2], 1) 
    alpha[1] <- 1 
    alpha[2] <- 1 
    alpha[3] <- 1 
    p.state[1:3] ~ ddirch(alpha[]) 

    for (t in 1:npts) { 
    idx[t] ~ dcat(p.state[]) 
    NSD[t] ~ dnorm(mu[idx[t]], tau[idx[t]])T(left[idx[t]], right[idx[t]]) 
    } 
} 

Усеченного расплывчатое нормальное распределение должно быть примерно эквивалентно равномерным распределение. Мы можем сравнить плотность ядер образцов от dnorm(0, 0.000001)T(0, 1) и dunif(0, 1). Здесь я использую откосы из R, но результат для WinBUGS должен быть похож:

library(R2jags) 
M <- ' 
model { 
    y_tnorm ~ dnorm(0, 0.000001)T(0, 1) 
    y_unif ~ dunif(0, 1) 
} 
' 
out <- jags(list(), NULL, c('y_tnorm', 'y_unif'), textConnection(M), 1, 100000, 
      n.burnin=0, n.thin=1, DIC=FALSE) 

plot(density(out$BUGSoutput$sims.matrix[, 'y_tnorm'], bw=0.001), main='') 
lines(density(out$BUGSoutput$sims.matrix[, 'y_unif'], bw=0.001), col=2) 
legend('bottomright', c('Truncated normal', 'Uniform'), bty='n', 
     col=1:2, lty=1, inset=0.05) 

---

EDIT

модели, кажется, прекрасно работать в Зубцове.

M <- 'model { 
    mu[1] ~ dunif(0, 1) 
    mu[2] <- min(mu[1] + eps[1], 1) 
    mu[3] <- 0.5 
    eps[1] ~ dunif(0, 1) 
    eps[2] ~ dunif(0, 1) 
    sigma[1] ~ dunif(0,.5)  
    sigma[2] ~ dunif(0,.5)  
    tau[1] <- pow(sigma[1],-2) 
    tau[2] <- pow(sigma[2],-2) 
    tau[3] <- 0.000001 
    left[1] <- -100 # something relatively very low 
    left[2] <- -100 # something relatively very low 
    left[3] ~ dunif(0, 1) 
    right[1] <- 100 # something relatively very high 
    right[2] <- 100 # something relatively very high 
    right[3] <- min(left[3] + eps[2], 1) 
    alpha[1] <- 1 
    alpha[2] <- 1 
    alpha[3] <- 1 
    p.state[1:3] ~ ddirch(alpha[]) 

    for (t in 1:npts) { 
    idx[t] ~ dcat(p.state[]) 
    NSD[t] ~ dnorm(mu[idx[t]], tau[idx[t]])T(left[idx[t]], right[idx[t]]) 
    } 
}' 


d <- read.csv('NSD.csv') 

library(R2jags) 
jagsfit <- jags(list(NSD=d$NSD, npts=nrow(d)), NULL, 
       c('mu', 'sigma', 'eps', 'left', 'right', 'p.state'), 
       textConnection(M), 3, 50000) 

Я не позволить ей работать достаточно долго, чтобы все параметры полностью сходятся, но вот предварительный взгляд на некоторые из ваших параметров.

##     mean  sd  2.5%  25%  50%  75%  97.5% Rhat n.eff 
## deviance -2650.2912 16.7002 -2667.7334 -2663.5577 -2656.7462 -2639.8387 -2610.2082 1.0054 450 
## eps[1]   0.9514 0.0021  0.9472  0.9500  0.9514  0.9528  0.9556 1.0018 2500 
## eps[2]   0.9100 0.0523  0.8438  0.8590  0.9018  0.9569  0.9975 1.0087 260 
## left[1]  -100.0000 0.0000 -100.0000 -100.0000 -100.0000 -100.0000 -100.0000 1.0000  1 
## left[2]  -100.0000 0.0000 -100.0000 -100.0000 -100.0000 -100.0000 -100.0000 1.0000  1 
## left[3]  0.0021 0.0013  0.0001  0.0011  0.0021  0.0032  0.0043 1.0011 14000 
## mu[1]   0.0008 0.0001  0.0007  0.0008  0.0008  0.0008  0.0009 1.0011 22000 
## mu[2]   0.9522 0.0021  0.9480  0.9508  0.9522  0.9536  0.9564 1.0017 2600 
## mu[3]   0.5000 0.0000  0.5000  0.5000  0.5000  0.5000  0.5000 1.0000  1 
## p.state[1]  0.4721 0.0259  0.4217  0.4546  0.4721  0.4898  0.5227 1.0010 60000 
## p.state[2]  0.3712 0.0265  0.3193  0.3532  0.3711  0.3890  0.4234 1.0017 2900 
## p.state[3]  0.1567 0.0207  0.1189  0.1423  0.1558  0.1700  0.1999 1.0019 2300 
## right[1]  100.0000 0.0000 100.0000 100.0000 100.0000 100.0000 100.0000 1.0000  1 
## right[2]  100.0000 0.0000 100.0000 100.0000 100.0000 100.0000 100.0000 1.0000  1 
## right[3]  0.9121 0.0522  0.8465  0.8610  0.9038  0.9589  0.9997 1.0087 260 
## sigma[1]  0.0007 0.0000  0.0006  0.0007  0.0007  0.0007  0.0008 1.0010 60000 
## sigma[2]  0.0238 0.0016  0.0210  0.0227  0.0238  0.0248  0.0272 1.0016 3200