Вы можете попробовать использовать неинформативный 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)
. Здесь я использую JAGS из 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)
![введите здесь описание изображения](https://i.stack.imgur.com/sjCfz.png)
ИЗМЕНИТЬ
Модель отлично работает в JAGS.
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
person
jbaums
schedule
19.02.2015
x[t, idx[t]]
. - person jbaums   schedule 19.02.2015