## =========================================================
## 4-component mixture-Normal prior
## Define the prior ONCE here
## =========================================================

m1 <- 0.0
m2 <- 0.0
m3 <- 0.0
m4 <- 0.0

sd1 <- .61
sd2 <- 1.42
sd3 <- 2.16
sd4 <- 5.64

w1 <- 0.32
w2 <- 0.31
w3 <- 0.3
w4 <- 0.07

## optional check: weights should sum to 1
sum(c(w1, w2, w3, w4))


## =========================================================
## p-values to analyse
## =========================================================

p_vals <- c(0.500, 0.300, 0.100, 0.050, 0.030, 0.010, 0.005)


## =========================================================
## Function to compute Pr(E > 0 | z)
## Model:
##   Z | E ~ N(E, 1)
## Prior:
##   E ~ sum_k w_k N(m_k, sd_k^2), k = 1,...,4
## =========================================================

mixnorm_postprob_table <- function(p_vals,
                                   m1, m2, m3, m4,
                                   sd1, sd2, sd3, sd4,
                                   w1, w2, w3, w4) {

  means   <- c(m1, m2, m3, m4)
  sds     <- c(sd1, sd2, sd3, sd4)
  weights <- c(w1, w2, w3, w4)

  weights <- weights / sum(weights)

  p_to_z <- function(p) {
    qnorm(1 - p / 2)
  }

  post_prob_gt0 <- function(z, means, sds, weights) {
    prior_var <- sds^2
    like_var  <- 1

    post_var  <- 1 / (1 / prior_var + 1 / like_var)
    post_mean <- post_var * (means / prior_var + z / like_var)

    marg_like <- dnorm(z, mean = means, sd = sqrt(prior_var + like_var))

    post_w_unnorm <- weights * marg_like
    post_w <- post_w_unnorm / sum(post_w_unnorm)

    comp_prob_gt0 <- 1 - pnorm(0, mean = post_mean, sd = sqrt(post_var))

    sum(post_w * comp_prob_gt0)
  }

  z_vals <- p_to_z(p_vals)

  pr_gt0 <- sapply(z_vals, function(z) {
    post_prob_gt0(z, means, sds, weights)
  })

  data.frame(
    p_value = p_vals,
    z_value = z_vals,
    Pr_E_gt_0_given_z = pr_gt0
  )
}


## =========================================================
## Run it
## =========================================================

results_table <- mixnorm_postprob_table(
  p_vals = p_vals,
  m1 = m1, m2 = m2, m3 = m3, m4 = m4,
  sd1 = sd1, sd2 = sd2, sd3 = sd3, sd4 = sd4,
  w1 = w1, w2 = w2, w3 = w3, w4 = w4
)

results_table