# initialization
mu.H0 <- 0; sigma.H0 <- 1
mu.H1 <- 0; sigma.H1 <- 1e-3
p.H1  <- 1/2
mu    <- c(mu.H0, mu.H1)
sigma <- c(sigma.H0, sigma.H1)

# simulation function
simulate <- function() {
  M <- rbinom(1, 1, p.H1); x <- rnorm(1, mu[M+1], sigma[M+1])
  x <- rnorm(1, mu[M+1], sigma[M+1])
#  p.val <- 2 * (1 - pnorm(mu[1] + abs(x-mu[1]), mu[1], sigma[1]) )
  p.val <- 2 * pnorm(mu[1] - abs(x-mu[1]), mu[1], sigma[1]) 
  BF <- dnorm(x, mu[2], sigma[2]) / dnorm(x, mu[1], sigma[1])
  lBF <- dnorm(x, mu[2], sigma[2], log=TRUE) - dnorm(x, mu[1], sigma[1], log=TRUE)
  cat(sprintf("x = %.5f  =>  p.val = %.2e,  BF = %.2e  [ log(BF) = %.2e ]\n",
              x, p.val, BF, lBF))
  return(M)
}