Bayesian Methods for the Physical Sciences by Stefano Andreon & Brian Weaver

Author:Stefano Andreon & Brian Weaver
Language: eng
Format: epub
Publisher: Springer International Publishing, Cham

7.7.3 Probabilistic Selection on the True Value

In this section, we would like to implement working code for the model illustrated in Sect. 7.5:

symbolic-model {

for (i in 1:length(obslnL)) {

obslnL[i] ~ dnorm(lnL[i],pow(err.lnL[i],-2)

lnL[i] ~ dnorm(lnLhat,pow(sigma.intr,-2)) phi((x-3)/2)

}

sigma.intr <- 3

lnLhat ~ dunif(-10,10)

}.

The product of a normal and a ϕ function is not among the JAGS distributions and it is not allowed in JAGS to write a line with a tilde and a product of a distribution and a function. The brute force approach is therefore compelling.

The current model is almost identical to the previous one, except for the fact that we replaced the previous sharp selection with the soft one. The latter replacement enters in three places: when we compute the soft selection p(lnL) (the correct[k] line), in the integral computation, and in the numerator of the posterior. For our comfort, we also changed the name of the variables. The model reads3:

data

{

# grid for integral evaluation

for (k in 1:400){

grid.x[k] <- (k-200.)/20.

correcc[k] <- phi((grid.x[k]-3)/2)

}

step.grid.x <-grid.x[2]-grid.x[1]

# dummy variable for zero-trick

# to sample from a distribution not available in JAGS

for (i in 1:length(obslnL)) {

dummy[i] <-0

}

}

model {

# compute the integral at denominator

for (k in 1:length(grid.x)) {

normaliz[k] <- sqrt(prec.intr/6.28)*exp(-prec.intr*

pow(lnLhat-grid.x[k],2 )/2.)*correcc[k]

}

tot.norm <- sum(normaliz[])*step.grid.x

for (i in 1:length(obslnL)) {

obslnL[i] ~ dnorm(lnL[i],prec.lnL[i])

# the following would suffice, if no selection

lnL[i] ~ dnorm(lnLhat,prec.intr)

# but selection is there. Adding the additional likelihood term

numerator[i] <- phi((lnL[i]-3)/2)

# sampling from an unavailable distribution

loglike[i] <- -log(numerator[i]/tot.norm)+C

dummy[i] ~ dpois(loglike[i])

}

C <- 10.

sigma.intr <- 3

prec.intr <- pow(sigma.intr,-2)

lnLhat ~ dunif(-10,10)

}.

Download

Copyright Disclaimer:
This site does not store any files on its server. We only index and link to content provided by other sites. Please contact the content providers to delete copyright contents if any and email us, we'll remove relevant links or contents immediately.