toscano2010_sim2.Rmd

---
title: "simu"
output:
  pdf_document: default
  html_document: default
date: "2026-05-26"
---

Run simulation 2 for Toscano & McMurray (2010).

Let $x$ be the observed cue value for a given cue dimension (e.g. VOT);

Let $i \in \{1, \dots, K\}$ index the K Gaussian components that together form the learner's internal category representation of a given cue dimension;

Let $\phi_i$ be the mixing weight, frequency of occurrence of the i-th Gaussian;

Let $\mu_i$ be its mean cue value of the Gaussian component along the cue dimension;

Let $\sigma_i$ be its standard deviation of the Gaussian component along the cue dimension;

Let $G_i(x)$ be the likelihood that the i-th Gaussian generated the observation $x$ given that cue dimension,

Let $\eta_\phi$, $\eta_\mu$, and $\eta_\sigma$ be the learning rates for the trial-by-trial updates of $\Delta\phi_i$, $\Delta\mu_i$, and $\Delta\sigma_i$ to the frequency, mean, and standard deviation of each Gaussian

$$
G_i(x) = \phi_i \frac{1}{\sqrt{2\pi\sigma_i^2}} \exp\left(-\frac{(x-\mu_i)^2}{2\sigma_i^2}\right) \tag{4}
$$

$$
M(x) = \sum_{i}^{K} G_i(x) \tag{5}
$$

Learning rule:

- A crucial innovation from this model is the use of a winner-take-all update rule such that only one Gaussian updates its value on any given trial. This allows the model to suppress unneeded categories and arrive at the correct solution. Without it, the model does not determine the correct number of categories (see McMurray et al., 2009a).

$$\Delta\phi_i = \eta_\phi \frac{G_i(x)}{M(x)} \tag{A2} $$

$$\Delta\mu_i = \eta_\mu \left(\frac{G_i(x)}{M(x)}\right) \frac{(x-\mu_i)}{\sigma_i^2} \tag{A3} $$

$$\Delta\sigma_i = \eta_\sigma \left(\frac{G_i(x)}{M(x)}\right) \left(\sigma_i^{-3}(x-\mu_i)^2 - \sigma_i^{-1}\right) \tag{A4} $$

# 1 generate data according to table 1

```{r}

library(ggplot2)
library(dplyr)
library(tidyr)

#synthesize data according to table 1
d.vot.voiced <- rnorm(n = 100000, mean = 0, sd = 5)
d.vl.voiced <- rnorm(n = 100000, mean = 188, sd = 45)
d.third.voiced <- rnorm(n = 100000, mean = 260, sd = 10)

d.vot.voiceless <- rnorm(n = 100000, mean = 50, sd = 10)
d.vl.voiceless <- rnorm(n = 100000, mean = 170, sd = 44)
d.third.voiceless <- rnorm(n = 100000, mean = 300, sd = 10)

d.voiced <- data.frame(category = "voiced",VOT= d.vot.voiced, 
                       VL = d.vl.voiced, Third  = d.third.voiced)

d.voiceless <- data.frame(category = "voiceless", VOT = d.vot.voiceless, 
                          VL = d.vl.voiceless, Third = d.third.voiceless)

d <- rbind(d.voiced, d.voiceless)

```

# 2 build up the model architecture for a single cue (e.g. VOT)

```{r}

mog.object <- function(K = 20, mu.m = 25, mu.sd = 75, sig = 3, seed = 42){
  list(mu = rnorm(K, mu.m, mu.sd),
    sigma = rep(sig, K),
    phi = rep(1 / K, K))
}

# Get G(x) and M(x) according to equation 4 and 5
mog.G <- function(mog.object, x) {
  mog.object$phi * dnorm(x, mean = mog.object$mu, sd = mog.object$sigma)
}
mog.M <- function(mog.object, x) sum(mog.G(mog.object, x))

mog.post <- function(mog.object, x){
  G <- mog.G(mog.object, x)
  M <- mog.M(mog.object, x)
  G/M
}

# learning rule according to A2-A5
# eta is for learning rate
# q.what is the use of sigma.min, sigma.max

mog.update <- function(mog.object, x,
                       eta.mu = 1, eta.sigma = 1, eta.phi = 0.001,
                       sigma.min = 0.1, sigma.max = NA){
  r <- mog.post(mog.object, x)
  mu.old <- mog.object$mu
  sg.old <- mog.object$sigma
  # A3 
  mog.object$mu <- mu.old + eta.mu * r * (x - mu.old) / sg.old^2
  # A4
  sig.new <- sg.old + eta.sigma * r * (sg.old^(-3) * (x - mu.old)^2 - sg.old^(-1))
  if (!is.na(sigma.max)) sig.new <- pmin(sig.new, sigma.max)
  mog.object$sigma <- pmax(sig.new, sigma.min)
  winner <- which.max(r)
  mog.object$phi[winner] <- mog.object$phi[winner] + eta.phi * r[winner]
  mog.object$phi <- pmax(mog.object$phi, 0)
  mog.object$phi <- mog.object$phi / sum(mog.object$phi)
  mog.object
}

# calculate weight, equation 6, using outer product
mog.weight <- function(mog.object){
  diffs <- outer(mog.object$mu, mog.object$mu, function(a, b) (a - b)^2)
  sigs  <- outer(mog.object$sigma, mog.object$sigma, "*")
  phis  <- outer(mog.object$phi, mog.object$phi, "*")
  0.5*sum(phis * diffs / sigs)
}

# sanity check on VOT only updating
mog <- mog.object(K = 20, mu.m = 25, mu.sd = 75, sig = 3)
idx <- sample(nrow(d), 90000, replace = TRUE) # shuffle so voiced/voiceless aren't blocked
for (i in idx) 
  mog <- mog.update(mog, d$VOT[i])

mog.weight(mog)

```

# 3 extend one cue to three cues (VOT, VL, combined)

p446. Initial $\mu$ values were randomly chosen from a distribution with a mean of 25 and standard deviation of 75 for the VOT dimension and a mean of 179 and standard deviation of 75 for the VL dimension. Initial $\sigma$ were set to 3 for the VOT dimiension and 10 for the VL dimension. K was set to 20, and initial $\phi$ were set to 1 ⁄ K. Learning rates were set to 1 ($\eta_\mu$), 1 ($\eta_\sigma$), 0.001 ($\eta_\phi$), and 0.001 ($\eta_p$).2 The models were then tested on a range of VOTs (0–40 ms in 5 ms steps) and two VLs (125 and 225 ms).

```{r}
vot.mog  <- mog.object(K = 20, mu.m = 25,  mu.sd = 75, sig = 3)
vot.mog.raw  <- mog.object(K = 20, mu.m = 25,  mu.sd = 75, sig = 3)

vl.mog <- mog.object(K = 20, mu.m = 179, mu.sd = 75, sig = 10)
comb.mog <- mog.object(K = 20, mu.m = 0,   mu.sd = 1,  sig = 1) 

vot.gm <- mean(d$VOT); vot.gs <- sd(d$VOT)
vl.gm  <- mean(d$VL);  vl.gs  <- sd(d$VL)
w.vot <- 0.5
w.vl  <- 0.5

idx <- sample(nrow(d), 5000, replace = TRUE)
for (i in idx) {
  vot.val <- d$VOT[i]
  vl.val  <- d$VL[i]
  # update the two cue MOGs from initial values
  vot.mog <- mog.update(vot.mog, vot.val,sigma.min = 0.1, sigma.max = 30)
  vl.mog  <- mog.update(vl.mog,  vl.val, sigma.min = 0.1, sigma.max = 80)
  # normalize the weights so that they sum up to one
  w.vot.raw <- mog.weight(vot.mog)
  w.vl.raw  <- mog.weight(vl.mog)
  tot <- w.vot.raw + w.vl.raw
  w.vot <- w.vot.raw / tot;  w.vl <- w.vl.raw / tot 
  # zscore the cue
  z.vot <- (vot.val - vot.gm) / vot.gs
  z.vl  <- (vl.val  - vl.gm ) / vl.gs
  x.comb <- w.vot * z.vot - w.vl * z.vl
  # update the combined ver as well, and pay attention to the value here, so that it won't explode
  comb.mog <- mog.update(comb.mog, w.vot * z.vot - w.vl * z.vl,
  eta.mu    = 0.01,
  eta.sigma = 0.01,
  eta.phi   = 0.001,
  sigma.min = 0.05,
  sigma.max = 1.0
)
}

```

# 4 resample and run multiple models

```{r}

# ---- setup (runs ONCE, before the loop) ----
n.runs   <- 100
n.trials <- 10000

mu.results      <- vector("list", n.runs)
comb.results    <- vector("list", n.runs)
weights.results <- data.frame(run = integer(), w.vot = numeric(), w.vl = numeric())
vot.results <- vector("list", n.runs)
vl.results  <- vector("list", n.runs)

posterior.p <- function(mog, x) {
  G  <- mog.G(mog, x)
  M  <- sum(G)
  gm <- sum(mog$phi * mog$mu)
  p.idx <- which.max(mog$phi * (mog$mu > gm))
  G[p.idx] / M
}

# ---- run loop ----
for (run in 1:n.runs) {
  set.seed(100 + run)

  vot.mog  <- mog.object(K = 20, mu.m = 25,  mu.sd = 75, sig = 3)
  vl.mog   <- mog.object(K = 20, mu.m = 179, mu.sd = 75, sig = 10)
  comb.mog <- mog.object(K = 20, mu.m = 0,   mu.sd = 1,  sig = 1)
  w.vot <- 0.1; w.vl <- 0.9

  idx <- sample(nrow(d), n.trials, replace = TRUE)
  for (i in idx) {
    vot.val <- d$VOT[i]; vl.val <- d$VL[i]
    vot.mog <- mog.update(vot.mog, vot.val)
    vl.mog  <- mog.update(vl.mog,  vl.val)
    w.vot.raw <- mog.weight(vot.mog); w.vl.raw <- mog.weight(vl.mog)
    tot <- w.vot.raw + w.vl.raw
    w.vot <- w.vot.raw / tot; w.vl <- w.vl.raw / tot
    z.vot <- (vot.val - vot.gm) / vot.gs
    z.vl  <- (vl.val  - vl.gm ) / vl.gs
    comb.mog <- mog.update(comb.mog, w.vot * z.vot - w.vl * z.vl, eta.mu    = 0.01,eta.sigma = 0.01,eta.phi   = 0.001,sigma.min = 0.05)
    #comb.mog <- mog.update(comb.mog, w.vot * z.vot - w.vl * z.vl)
  }

  vot.a  <- which(vot.mog$phi  > 0.1)
  vl.a   <- which(vl.mog$phi   > 0.1)
  comb.a <- which(comb.mog$phi > 0.1)
  converged <- (length(vot.a) == 2 &&length(vl.a)  == 2 &&length(comb.a) == 2)

  #cat(sprintf("Run %2d: VOT=%d VL=%d COMB=%d  %s\n",run, length(vot.a), length(vl.a), length(comb.a),ifelse(converged, "OK", "skip")))

  if (converged) {
    vot.ord <- order(vot.mog$mu[vot.a])                  # voiced first (smaller VOT)
    vl.ord  <- order(vl.mog$mu[vl.a], decreasing = TRUE) # voiced first (larger VL)

    mu.results[[run]] <- data.frame(
      run = run,
      cue = c("VOT","VOT","VL","VL"),
      category = c("voiced","voiceless","voiced","voiceless"),
      mu    = c(vot.mog$mu[vot.a][vot.ord],
                vl.mog$mu[vl.a][vl.ord]),
      sigma = c(vot.mog$sigma[vot.a][vot.ord],
                vl.mog$sigma[vl.a][vl.ord])
    )
    weights.results <- rbind(weights.results,
                             data.frame(run = run, w.vot = w.vot, w.vl = w.vl))
    vot.results[[run]] <- vot.mog
    vl.results[[run]]  <- vl.mog
    comb.results[[run]] <- comb.mog
  }
}  


params <- do.call(rbind, mu.results)
#cat(sprintf("\n%d / %d runs survived (>=2 categories in both cue MOGs)\n",length(unique(params$run)), n.runs))



```





# 5 Check the results


## Visualize trading relation

```{r}

avg.mog <- function(mog.list, threshold = 0.1) {
  mogs <- mog.list[!sapply(mog.list, is.null)]
  surv <- lapply(mogs, function(m) {
    a <- which(m$phi > threshold)
    o <- order(m$mu[a])                       # ascending mu
    list(phi = m$phi[a][o], mu = m$mu[a][o], sigma = m$sigma[a][o])
  })
  surv <- surv[sapply(surv, function(s) length(s$mu) == 2)]
  stopifnot(length(surv) > 0)
  list(
    phi   = colMeans(do.call(rbind, lapply(surv, `[[`, "phi"))),
    mu    = colMeans(do.call(rbind, lapply(surv, `[[`, "mu"))),
    sigma = colMeans(do.call(rbind, lapply(surv, `[[`, "sigma")))
  )
}

vot.avg  <- avg.mog(vot.results)
vl.avg   <- avg.mog(vl.results)
comb.avg <- avg.mog(comb.results)

avg.w.vot <- mean(weights.results$w.vot)
avg.w.vl  <- mean(weights.results$w.vl)

N <- 200
vot.r <- seq(-20,  80, length.out = N)
vl.r  <- seq( 50, 350, length.out = N)
cmb.r <- seq( -3,   3, length.out = N)
xn    <- seq(  0,   1, length.out = N)

post.df <- rbind(
  data.frame(xn = xn, p = sapply(vot.r, function(x) posterior.p(vot.avg,  x)), MOG = "VOT"),
  data.frame(xn = xn, p = sapply(vl.r,  function(x) posterior.p(vl.avg,   x)), MOG = "VL"),
  data.frame(xn = xn, p = sapply(cmb.r, function(x) posterior.p(comb.avg, x)), MOG = "Combined")
)
post.df$MOG <- factor(post.df$MOG, levels = c("VOT","VL","Combined"))

ggplot(post.df, aes(xn, p, colour = MOG)) +
  geom_line(linewidth = 1.3) +
  scale_colour_manual(values = c(VOT = "#7BB8E8", VL = "#FFB347", Combined = "#90C695")) +
  ylim(0, 1) +
  labs(x = "Cue value (normalized)", y = "P(/p/-side | x)",
       title = "Simulation 2 — Posteriors per MOG",
       subtitle = sprintf("Averaged over %d successful runs", nrow(weights.results))) +
  theme_minimal(base_size = 13) +
  theme(legend.position = "top")

vot.test <- seq(0, 40, by = 5)
fig7 <- expand.grid(VOT = vot.test, VL = c(125, 225))
fig7$p_p <- mapply(function(vot, vl) {
  z.vot <- (vot - vot.gm) / vot.gs
  z.vl  <- (vl  - vl.gm ) / vl.gs
  x.c   <- avg.w.vot * z.vot - avg.w.vl * z.vl
  posterior.p(comb.avg, x.c)
}, fig7$VOT, fig7$VL)
fig7$VL <- factor(fig7$VL, levels = c(125, 225),
                  labels = c("VL = 125 ms", "VL = 225 ms"))

ggplot(fig7, aes(VOT, p_p, colour = VL)) +
  geom_line(linewidth = 1.4) +
  geom_point(size = 2.5) +
  scale_colour_manual(values = c("VL = 125 ms" = "#7BB8E8",
                                 "VL = 225 ms" = "#FFB347")) +
  ylim(0, 1) +
  labs(x = "VOT (ms)", y = "Proportion /p/", title = "Simulation 2 — Trading relation (Fig. 7)", colour = NULL) +
  theme_minimal(base_size = 13) +
  theme(legend.position = "top")
```



## Compare learned and veridical parameters

```{r}

# ---- Table 1 reference: true mu, sigma, and Eq. 6 weights ----
true_table1 <- data.frame(
  cue      = c("VOT","VOT","VL","VL"),
  category = c("voiced","voiceless","voiced","voiceless"),
  mu_true    = c(0,  50, 188, 170),
  sigma_true = c(5,  10,  45,  44)
)

# Eq. 6 evaluated on Table 1 directly (two-category case, phi=0.5 each):
# w = phi_v * phi_u * (mu_v - mu_u)^2 / (sigma_v * sigma_u) = 0.25 * (.)^2 / (..)
w_t1_vot <- 0.25 * (0 - 50)^2  / (5  * 10)        # = 12.5
w_t1_vl  <- 0.25 * (188 - 170)^2 / (45 * 44)      # = 0.041
tot_t1   <- w_t1_vot + w_t1_vl
w_table1 <- data.frame(cue = c("VOT","VL"),
                        weight_true = c(w_t1_vot / tot_t1,
                                        w_t1_vl  / tot_t1))

# ---- Summarize learned mu & sigma across runs ----
summary_df <- params %>%
  filter(!is.na(mu))%>%
  group_by(cue, category) %>%
  summarise(mu.mean    = mean(mu),    mu.sd    = sd(mu),
            sigma.mean = mean(sigma), sigma.sd = sd(sigma),
            .groups = "drop") %>%
  left_join(true_table1, by = c("cue", "category"))


# ---- Plot 1: learned mu vs Table 1 mu ----
ggplot(summary_df, aes(x = category)) +
  geom_pointrange(aes(y = mu.mean, ymin = mu.mean - mu.sd, ymax = mu.mean + mu.sd,
                       colour = "Learned (mean ± SD across runs)"),
                   size = 0.7) +
  geom_point(aes(y = mu_true, colour = "Table 1"),
              shape = 4, size = 4, stroke = 1.5) +
  facet_wrap(~ cue, scales = "free_y") +
  scale_colour_manual(values = c("Learned (mean ± SD across runs)" = "steelblue",
                                  "Table 1" = "black")) +
  labs(title = "Learned vs Table 1 category means (μ)",
       y = "μ", x = NULL, colour = NULL) +
  theme_minimal()


# ---- Plot 2: learned sigma vs Table 1 sigma ----
ggplot(summary_df, aes(x = category)) +
  geom_pointrange(aes(y = sigma.mean, ymin = sigma.mean - sigma.sd, ymax = sigma.mean + sigma.sd,
                       colour = "Learned (mean ± SD across runs)"),
                   size = 0.7) +
  geom_point(aes(y = sigma_true, colour = "Table 1"),
              shape = 4, size = 4, stroke = 1.5) +
  facet_wrap(~ cue, scales = "free_y") +
  scale_colour_manual(values = c("Learned (mean ± SD across runs)" = "firebrick",
                                  "Table 1" = "black")) +
  labs(title = "Learned vs Table 1 category SDs (σ)",
       y = "σ", x = NULL, colour = NULL) +
  theme_minimal()


# ---- Plot 3: cue weights, one dot per run, vs Eq. 6 on Table 1 ----
w_long <- weights.results %>%
  pivot_longer(c(w.vot, w.vl), names_to = "cue", values_to = "weight") %>%
  mutate(cue = recode(cue, w.vot = "VOT", w.vl = "VL"))

w_summary <- w_long %>%
  group_by(cue) %>%
  summarise(w.mean = mean(weight), w.sd = sd(weight), .groups = "drop")

ggplot() +
  geom_jitter(data = w_long, aes(x = cue, y = weight, colour = "Learned (per run)"),
               width = 0.08, alpha = 0.45) +
  geom_pointrange(data = w_summary,
                   aes(x = cue, y = w.mean,
                       ymin = w.mean - w.sd, ymax = w.mean + w.sd,
                       colour = "Learned (mean ± SD)"),
                   size = 0.8) +
  geom_point(data = w_table1,
              aes(x = cue, y = weight_true, colour = "Table 1 (Eq. 6 on truth)"),
              shape = 4, size = 4, stroke = 1.5) +
  scale_colour_manual(values = c("Learned (per run)"      = "steelblue",
                                  "Learned (mean ± SD)"   = "steelblue4",
                                  "Table 1 (Eq. 6 on truth)" = "black")) +
  ylim(0, 1) +
  labs(title = "Cue weights: learned across runs vs Eq. 6 on Table 1",
       y = "normalized weight", x = NULL, colour = NULL) +
  theme_minimal()
```
Observations: 1) the unsupervised learning algorithm approximates the veridical mean and SD for both voiced and voiceless categories relatively well along VOT dimension, 2) but it has a tendency to overestimate mean and underestimate SD for both voiced and voiceless categories along VL dimension (-why), 3) the learned weight for VL and VOT also shrinks from the veridical weight. 



# 7 Compare against an ideal Bayesian model cue integration

Let $v$ = VOT observation,

Let $\ell$ = VL observation,

Let C ∈ {b, p} the category.


$$\mathcal N(x; \mu, \sigma^2) = \frac{1}{\sigma\sqrt{2\pi}} \exp\!\Big(-\tfrac{(x-\mu)^2}{2\sigma^2}\Big)$$

The likelihood for /p/ will be: 
$$L_p \;=\; P(v, \ell \mid C=p) \;=\; \mathcal N(v; \mu_{v,p}, \sigma_{v,p}^2)\;\cdot\;\mathcal N(\ell; \mu_{\ell,p}, \sigma_{\ell,p}^2)$$

The **Bayes' rule** for the posterior of /p/:

$$P(C=p \mid v, \ell) \;=\; \frac{P(v, \ell \mid p)\, P(p)}{P(v, \ell \mid p)\, P(p) \;+\; P(v, \ell \mid b)\, P(b)}$$
The denominator is the marginal likelihood (law of total probability) — it sums over the two possible categories so the posterior is normalized to 1. With flat prior P(p) = P(b) = 0.5, the prior factors cancel and you get the **likelihood ratio form**:

$$P(p \mid v, \ell) \;=\; \frac{L_p}{L_p + L_b}$$


Use the MAP decision rule

Choose /p/ iff P(p | v, ℓ) > 0.5. Equivalently, iff the **log posterior ratio** is positive:

$$\log \frac{P(p \mid v, \ell)}{P(b \mid v, \ell)} \;>\; 0$$

Expanding under conditional independence and flat prior, this decomposes into a per-cue sum:

$${\log\frac{\mathcal N(v;\mu_{v,p},\sigma_{v,p}^2)}{\mathcal N(v;\mu_{v,b},\sigma_{v,b}^2)}}_{\text{}} \;+\; {\log\frac{\mathcal N(\ell;\mu_{\ell,p},\sigma_{\ell,p}^2)}{\mathcal N(\ell;\mu_{\ell,b},\sigma_{\ell,b}^2)}}_{} \;>\; 0$$

Each log-likelihood ratio for a Gaussian works out to:

$$\log \frac{\mathcal N(x; \mu_p, \sigma_p^2)}{\mathcal N(x; \mu_b, \sigma_b^2)} \;=\; \log\frac{\sigma_b}{\sigma_p} \;-\; \frac{(x-\mu_p)^2}{2\sigma_p^2} \;+\; \frac{(x-\mu_b)^2}{2\sigma_b^2}$$

So the boundary in (v, ℓ) space is the level set of a sum of two quadratics. With equal within-cue variances per category the quadratic terms cancel and you get a straight line; with VOT's unequal variances (5 vs 10) the boundary is slightly curved.


These lines step outside the posterior calculation and compute the optimal cue weights under a simplifying equal-variance approximation.

 d-prime, the SDT discriminability index:

$$d'_i \;=\; \frac{\mu_{i,p} - \mu_{i,b}}{\sigma_{\text{pool},i}}$$

This measures category separation in units of within-category SD — i.e., a signal-to-noise ratio for binary discrimination. Plugging in:

$$d'_\text{VOT} = \frac{50 - 0}{7.91} \approx 6.32 \qquad d'_\text{VL} = \frac{170 - 188}{44.5} \approx -0.40$$
**`J.vot` and `J.vl`** — precision (inverse variance in d-prime units):

$$J_i \;=\; (d'_i)^2$$

For Gaussian likelihoods, d'² is proportional to the **Fisher information** the cue provides about the binary category. Equivalently, it's the inverse variance of the optimal unbiased category-indicator estimate from that cue (the y_i from the earlier breakdown).

$$J_\text{VOT} \approx 39.97 \qquad J_\text{VL} \approx 0.16$$

**`w.vot.opt` and `w.vl.opt`** — reliability weights, the precision-weighted-average formula:

$$w_i \;=\; \frac{J_i}{J_\text{VOT} + J_\text{VL}} \;=\; \frac{(d'_i)^2}{\sum_j (d'_j)^2}$$

These are the weights the textbook formula μ_post = (J_1 x_1 + J_2 x_2)/(J_1 + J_2) assigns to each cue when you recast the problem as estimating a common latent category indicator from two noisy reads. Plugging in:

$$w_\text{VOT} \approx \frac{39.97}{39.97 + 0.16} \approx 0.9960$$

$$w_\text{VL} \approx \frac{0.16}{39.97 + 0.16} \approx 0.0040$$

So the precision-weighted-average view says: an optimal Bayesian observer assigns ~99.6% of the cue weight to VOT and ~0.4% to VL — confirming the Fig. 8 baseline.

The `posterior.p.bayes` function does the **full** Bayesian computation (exact, including unequal variances and the curvature of the decision boundary). The d'-based weights are a **summary statistic** that compactly describes what the posterior is doing on average: under the equal-variance simplification, the log-posterior ratio is linear in each cue, and its slope along cue i is exactly (μ_p,i − μ_b,i)/σ²_i — which after rescaling each cue gives you the J_i / (J_1+J_2) weights.

So both blocks describe the same observer; the posterior function is what you'd actually compute for predictions, and the d'² weights are the diagnostic that tells you *why* the trading relation is absent — VL's d' is just too small.


```{r baseline}
# ---- Bayesian cue integration with veridical Table-1 parameters ----

# Category-conditional cue parameters (mean, sd)
mu.vot.b  <- 0;   sd.vot.b  <- 5
mu.vot.p  <- 50;  sd.vot.p  <- 10
mu.vl.b   <- 188; sd.vl.b   <- 45
mu.vl.p   <- 170; sd.vl.p   <- 44

prior.b <- 0.5
prior.p <- 0.5

# Analytic posterior P(/p/ | VOT, VL) under conditional independence,
# with exact (unequal-variance) Gaussian likelihoods per cue per category.
posterior.p.bayes <- function(vot, vl, prior.p = 0.5, prior.b = 0.5) {
  L.p <- dnorm(vot, mu.vot.p, sd.vot.p) * dnorm(vl, mu.vl.p, sd.vl.p)
  L.b <- dnorm(vot, mu.vot.b, sd.vot.b) * dnorm(vl, mu.vl.b, sd.vl.b)
  (L.p * prior.p) / (L.p * prior.p + L.b * prior.b)
}

# MAP decision rule (binary)
decide.bayes <- function(vot, vl, threshold = 0.5)
  posterior.p.bayes(vot, vl) > threshold

# ---- Veridical "optimal" cue weights from d-prime^2 (for reference) ----
sd.vot.pool <- sqrt((sd.vot.b^2 + sd.vot.p^2) / 2)
sd.vl.pool  <- sqrt((sd.vl.b^2  + sd.vl.p^2)  / 2)
dp.vot <- (mu.vot.p - mu.vot.b) / sd.vot.pool
dp.vl  <- (mu.vl.p  - mu.vl.b ) / sd.vl.pool
J.vot  <- dp.vot^2
J.vl   <- dp.vl^2
w.vot.opt <- J.vot / (J.vot + J.vl)
w.vl.opt  <- J.vl  / (J.vot + J.vl)
cat(sprintf("d' VOT=%.2f  d' VL=%.2f\n", dp.vot, dp.vl))
cat(sprintf("Optimal weights:  w_VOT=%.4f  w_VL=%.4f\n", w.vot.opt, w.vl.opt))
# Expect approx 0.996 / 0.004 — matches Toscano Fig. 8 (0.997/0.003)

# ---- Trading-relation grid (same test points as Fig. 7) ----
vot.test <- seq(0, 40, by = 5)
bayes.df <- expand.grid(VOT = vot.test, VL = c(125, 225))
bayes.df$p_p <- posterior.p.bayes(bayes.df$VOT, bayes.df$VL)
bayes.df$VL  <- factor(bayes.df$VL, levels = c(125, 225),
                       labels = c("VL = 125 ms", "VL = 225 ms"))

# ---- Stand-alone plot (Bayesian only) ----
ggplot(bayes.df, aes(VOT, p_p, colour = VL)) +
  geom_line(linewidth = 1.4) + geom_point(size = 2.5) +
  scale_colour_manual(values = c("VL = 125 ms" = "#7BB8E8",
                                 "VL = 225 ms" = "#FFB347")) +
  ylim(0, 1) +
  labs(x = "VOT (ms)", y = "P(/p/ | VOT, VL)",
       title = "Optimal Bayes — veridical parameters",
       subtitle = "Expect a near-zero trading relation (the Fig. 8 baseline)") +
  theme_minimal(base_size = 13) +
  theme(legend.position = "top")
```


```{r with d1 d2 }
# ---- Bayesian cue integration with veridical Table-1 parameters ----

# Category-conditional cue parameters (mean, sd)

d1=60
d2 =1

mu.vot.b  <- 0;   sd.vot.b  <- 5
mu.vot.p  <- 50;  sd.vot.p  <- 10
mu.vl.p   <- 170; sd.vl.p   <- 44
mu.vl.b <- mu.vl.p+ d1; sd.vl.p <- sd.vl.p +d2
#mu.vl.b <- 188; mu.vl.p <- 45


prior.b <- 0.5
prior.p <- 0.5

# Analytic posterior P(/p/ | VOT, VL) under conditional independence,
# with exact (unequal-variance) Gaussian likelihoods per cue per category.
posterior.p.bayes <- function(vot, vl, prior.p = 0.5, prior.b = 0.5) {
  L.p <- dnorm(vot, mu.vot.p, sd.vot.p) * dnorm(vl, mu.vl.p, sd.vl.p)
  L.b <- dnorm(vot, mu.vot.b, sd.vot.b) * dnorm(vl, mu.vl.b, sd.vl.b)
  (L.p * prior.p) / (L.p * prior.p + L.b * prior.b)
}

# MAP decision rule (binary)
decide.bayes <- function(vot, vl, threshold = 0.5)
  posterior.p.bayes(vot, vl) > threshold

# ---- Veridical "optimal" cue weights from d-prime^2 (for reference) ----
sd.vot.pool <- sqrt((sd.vot.b^2 + sd.vot.p^2) / 2)
sd.vl.pool  <- sqrt((sd.vl.b^2  + sd.vl.p^2)  / 2)
dp.vot <- (mu.vot.p - mu.vot.b) / sd.vot.pool
dp.vl  <- (mu.vl.p  - mu.vl.b ) / sd.vl.pool
J.vot  <- dp.vot^2
J.vl   <- dp.vl^2
w.vot.opt <- J.vot / (J.vot + J.vl)
w.vl.opt  <- J.vl  / (J.vot + J.vl)
cat(sprintf("d' VOT=%.2f  d' VL=%.2f\n", dp.vot, dp.vl))
cat(sprintf("Optimal weights:  w_VOT=%.4f  w_VL=%.4f\n", w.vot.opt, w.vl.opt))
# Expect approx 0.996 / 0.004 — matches Toscano Fig. 8 (0.997/0.003)

# ---- Trading-relation grid (same test points as Fig. 7) ----
vot.test <- seq(0, 40, by = 5)
bayes.df <- expand.grid(VOT = vot.test, VL = c(125, 225))
bayes.df$p_p <- posterior.p.bayes(bayes.df$VOT, bayes.df$VL)
bayes.df$VL  <- factor(bayes.df$VL, levels = c(125, 225),
                       labels = c("VL = 125 ms", "VL = 225 ms"))

# ---- Stand-alone plot (Bayesian only) ----
ggplot(bayes.df, aes(VOT, p_p, colour = VL)) +
  geom_line(linewidth = 1.4) + geom_point(size = 2.5) +
  scale_colour_manual(values = c("VL = 125 ms" = "#7BB8E8",
                                 "VL = 225 ms" = "#FFB347")) +
  ylim(0, 1) +
  labs(x = "VOT (ms)", y = "P(/p/ | VOT, VL)",
       title = "Optimal Bayes — veridical parameters",
       subtitle = "Expect a near-zero trading relation (the Fig. 8 baseline)") +
  theme_minimal(base_size = 13) +
  theme(legend.position = "top")
```
Make the following:

d1 = VL category mean separation (μ_VL,b − μ_VL,p)

d2 = VL within-category SD (shared across /b/ and /p/ for the sweep)

```{r}

# ---- Fixed parameters ----
mu.vot.b <- 0;   sd.vot.b <- 5
mu.vot.p <- 50;  sd.vot.p <- 10
mu.vl.p  <- 170                # /p/ VL mean fixed; d1 sets /b/ separation, d2 sets VL SD

# ---- One (d1, d2) -> trading-relation curve ----
trading.curve <- function(d1, d2,
                          vot.test = seq(0, 40, by = 5),
                          vl.test  = c(125, 225)) {
  mu.vl.b.d <- mu.vl.p + d1
  sd.vl     <- d2

  post <- function(vot, vl) {
    L.p <- dnorm(vot, mu.vot.p, sd.vot.p) * dnorm(vl, mu.vl.p,    sd.vl)
    L.b <- dnorm(vot, mu.vot.b, sd.vot.b) * dnorm(vl, mu.vl.b.d, sd.vl)
    L.p / (L.p + L.b)
  }

  df <- expand.grid(VOT = vot.test, VL = vl.test)
  df$p_p       <- post(df$VOT, df$VL)
  df$d1        <- d1
  df$d2        <- d2
  df$delta.llr <- (-d1) / d2^2 * (125 - 225)   # equal-variance closed form
  df
}

# ---- Sweep over a grid of (d1, d2) ----
d1.values <- c(5, 18, 40, 80)      # VL separation: small → veridical → exaggerated → huge
d2.values <- c(10, 25, 44, 80)     # VL SD: tight → medium → veridical → wide

grid       <- expand.grid(d1 = d1.values, d2 = d2.values)
all.curves <- do.call(rbind, Map(trading.curve, grid$d1, grid$d2))

all.curves <- all.curves %>%
  mutate(
    VL.label = factor(VL, levels = c(125, 225),
                      labels = c("VL = 125 ms", "VL = 225 ms")),
    d1.label = factor(sprintf("d1 = %d ms", d1),
                      levels = sprintf("d1 = %d ms", d1.values)),
    d2.label = factor(sprintf("d2 = %d ms", d2),
                      levels = sprintf("d2 = %d ms", d2.values))
  )

# ---- Trading-relation panel grid ----
ggplot(all.curves, aes(VOT, p_p, colour = VL.label)) +
  geom_line(linewidth = 1) +
  geom_point(size = 1.4) +
  facet_grid(d2.label ~ d1.label) +
  scale_colour_manual(values = c("VL = 125 ms" = "#7BB8E8",
                                 "VL = 225 ms" = "#FFB347")) +
  ylim(0, 1) +
  labs(x = "VOT (ms)", y = "P(/p/ | VOT, VL)",
       title = "Trading relation vs. VL category separation (d1) and VL SD (d2)",
       subtitle = "Rows = d2 (within-category SD), Columns = d1 (category mean separation)",
       colour = NULL) +
  theme_minimal(base_size = 11) +
  theme(legend.position = "top")
```
```{r parameter estimation}

# ---- synthetic listener target (replace with real data when you have it) ----
# Pretend listeners produce a trading relation with Δ_LLR ≈ 2 log-odds (≈ 0.5 gap)
listener.target <- function(vot, vl, k_VL = 22) {
  # generated from the model with exaggerated VL slope; used as stand-in for real data
  sd <- 44; mu.b <- 170 + k_VL; mu.p <- 170
  L.p <- dnorm(vot, 50, 10) * dnorm(vl, mu.p, sd)
  L.b <- dnorm(vot, 0,  5)  * dnorm(vl, mu.b, sd)
  L.p / (L.p + L.b)
}

vot.grid <- seq(0, 40, by = 5)
vl.grid  <- c(125, 225)
test     <- expand.grid(VOT = vot.grid, VL = vl.grid)
n.per.cell <- 50                                # trials per cell (set to your N)
test$p_listener <- listener.target(test$VOT, test$VL)
test$y          <- rbinom(nrow(test), size = n.per.cell, prob = test$p_listener)

# ---- model prediction at given (d1, d2) ----
model.p <- function(vot, vl, d1, d2) {
  mu.vl.b.d <- 170 + d1
  L.p <- dnorm(vot, 50, 10) * dnorm(vl, 170,        d2)
  L.b <- dnorm(vot, 0,  5)  * dnorm(vl, mu.vl.b.d,  d2)
  L.p / (L.p + L.b)
}

# ---- negative log-likelihood ----
nll <- function(par) {
  d1 <- par[1]; d2 <- par[2]
  if (d2 <= 0) return(1e10)
  p <- model.p(test$VOT, test$VL, d1, d2)
  p <- pmin(pmax(p, 1e-9), 1 - 1e-9)
  -sum(dbinom(test$y, n.per.cell, p, log = TRUE))
}

fit <- optim(par = c(d1 = 20, d2 = 40), fn = nll,
             method = "L-BFGS-B",
             lower = c(0.1, 1), upper = c(200, 200))

cat(sprintf("MLE:  d1 = %.2f,  d2 = %.2f\n", fit$par[1], fit$par[2]))
cat(sprintf("VL slope (d1/d2^2) = %.4f log-odds/ms\n", fit$par[1] / fit$par[2]^2))

```
```{r}

# (assumes you've already run Framing A and have: test, n.per.cell, fit, model.p)

# ---- 1. Fitted curves vs. listener data ----
vot.fine <- seq(0, 40, length.out = 100)
fit.df <- expand.grid(VOT = vot.fine, VL = vl.grid)
fit.df$p_p <- model.p(fit.df$VOT, fit.df$VL, fit$par[1], fit$par[2])

# binomial confidence intervals for the data

ci <- function(y, n) {
  out <- t(sapply(seq_along(y), function(i) {
    bt <- binom.test(y[i], n)
    c(p = unname(bt$estimate), lo = bt$conf.int[1], hi = bt$conf.int[2])
  }))
  as.data.frame(out)
}

data.df      <- cbind(test, ci(test$y, n.per.cell))
data.df.both <- data.df

fit.df$VL.label  <- factor(fit.df$VL,  labels = c("VL = 125 ms","VL = 225 ms"))
data.df$VL.label <- factor(data.df$VL, labels = c("VL = 125 ms","VL = 225 ms"))

p.fit <- ggplot() +
  geom_ribbon(data = fit.df, aes(VOT, ymin = NA, ymax = NA), inherit.aes = FALSE) +
  geom_line(data = fit.df, aes(VOT, p_p, colour = VL.label), linewidth = 1.2) +
  geom_pointrange(data = data.df,
                  aes(VOT, p, ymin = lo, ymax = hi, colour = VL.label),
                  size = 0.4) +
  scale_colour_manual(values = c("VL = 125 ms" = "#7BB8E8",
                                 "VL = 225 ms" = "#FFB347")) +
  ylim(0, 1) +
  labs(x = "VOT (ms)", y = "Proportion /p/",
       title = "Framing A — MLE fit to listener data",
       subtitle = sprintf("MLE: d1 = %.1f, d2 = %.1f  (VL slope = %.4f)",
                          fit$par[1], fit$par[2],
                          fit$par[1] / fit$par[2]^2),
       colour = NULL) +
  theme_minimal(base_size = 13) +
  theme(legend.position = "top")
print(p.fit)

# ---- 2. Negative log-likelihood surface (identifiability ridge) ----
d1.range <- seq(5, 80, length.out = 60)
d2.range <- seq(10, 80, length.out = 60)
surf <- expand.grid(d1 = d1.range, d2 = d2.range)
surf$nll <- mapply(function(a, b) nll(c(a, b)), surf$d1, surf$d2)
surf$rel.nll <- surf$nll - min(surf$nll)        # relative to MLE

p.surf <- ggplot(surf, aes(d1, d2, fill = pmin(rel.nll, 50))) +
  geom_raster() +
  geom_contour(aes(z = pmin(rel.nll, 50)),
               colour = "white", linewidth = 0.3,
               breaks = c(2, 5, 10, 20, 40)) +
  # mark the MLE
  annotate("point", x = fit$par[1], y = fit$par[2],
           colour = "red", size = 4, shape = 4, stroke = 1.5) +
  annotate("text",  x = fit$par[1], y = fit$par[2] + 4,
           label = "MLE", colour = "red", hjust = 0.5) +
  # mark the veridical
  annotate("point", x = 18, y = 44, colour = "white", size = 3) +
  annotate("text",  x = 22, y = 44, label = "veridical",
           colour = "white", hjust = 0) +
  # overlay the d1 = c * d2^2 ridge (constant VL slope = MLE slope)
  stat_function(fun = function(d2) (fit$par[1] / fit$par[2]^2) * d2^2,
                colour = "yellow", linetype = "dashed", linewidth = 0.8,
                inherit.aes = FALSE) +
  scale_fill_viridis_c(name = "ΔNLL", option = "magma", direction = -1) +
  coord_cartesian(xlim = range(d1.range), ylim = range(d2.range)) +
  labs(x = "d1 (VL separation, ms)", y = "d2 (VL SD, ms)",
       title = "Likelihood surface: the d1/d2² identifiability ridge",
       subtitle = "Yellow dashed = constant VL slope; X = MLE; • = veridical (acoustic) values") +
  theme_minimal(base_size = 12)
print(p.surf)

# ---- 3. Veridical vs. MLE side-by-side ----
veridical.p <- model.p(fit.df$VOT, fit.df$VL, d1 = 18, d2 = 44)
compare.df <- rbind(
  data.frame(fit.df, model = "MLE fit"),
  data.frame(VOT = fit.df$VOT, VL = fit.df$VL, VL.label = fit.df$VL.label,
             p_p = veridical.p, model = "Veridical Bayes")
)
data.df.both <- data.df

p.compare <- ggplot() +
  geom_line(data = compare.df,
            aes(VOT, p_p, colour = VL.label, linetype = model),
            linewidth = 1.1) +
  geom_pointrange(data = data.df.both,
                  aes(VOT, p.p, ymin = lo, ymax = hi, colour = VL.label),
                  size = 0.35, alpha = 0.8) +
  scale_colour_manual(values = c("VL = 125 ms" = "#7BB8E8",
                                 "VL = 225 ms" = "#FFB347")) +
  scale_linetype_manual(values = c("MLE fit" = "solid",
                                   "Veridical Bayes" = "dashed")) +
  ylim(0, 1) +
  labs(x = "VOT (ms)", y = "Proportion /p/",
       title = "Listener data vs. veridical Bayes vs. MLE-fit Bayes",
       colour = NULL, linetype = NULL) +
  theme_minimal(base_size = 13) +
  theme(legend.position = "top", legend.box = "vertical")
print(p.compare)
```


# appendix: the derivation of A2

We want $\frac{\partial \log M(x)}{\partial \mu_i}$. 

Apply the log rule

$$
\frac{\partial \log M}{\partial \mu_i} = \frac{1}{M} \cdot \frac{\partial M}{\partial \mu_i}  = \frac{1}{M} \cdot \sum_j \frac{\partial G_j}{\partial \mu_i} 
$$

And as $G_j$ only contains $\mu_j$, not the other $\mu$'s. So $\frac{\partial G_j}{\partial \mu_i} = 0$ and only the $j=i$ term survives

$$
\frac{\partial G_j}{\partial \mu_i} = \frac{\partial G_i}{\partial \mu_i} = \phi_i \cdot \frac{1}{\sqrt{2\pi\sigma_i^2}} \cdot \frac{\partial}{\partial \mu_i}\left[ \exp!\left( -\frac{(x-\mu_i)^2}{2\sigma_i^2} \right) \right]
$$

Apply the exponential rule and the chain rule we get the following

$$
u(\mu_i) = -\frac{(x-\mu_i)^2}{2\sigma_i^2} \\[3pt]\frac{du}{d\mu_i} = -\frac{1}{2\sigma_i^2} \cdot \frac{d}{d\mu_i}\left[(x-\mu_i)^2\right]
$$

apply chain rule once again

$$
\frac{d}{d\mu_i}(x - \mu_i)^2 = 2(x - \mu_i) \cdot \frac{d}{d\mu_i}(x - \mu_i) = 2(x - \mu_i) \cdot (-1) = -2(x - \mu_i) \\[3pt]
\frac{du}{d\mu_i} = -\frac{1}{2\sigma_i^2} \cdot (-2)(x - \mu_i) = \frac{x - \mu_i}{\sigma_i^2}
$$

and plug it back
$$
\frac{\partial G_i}{\partial \mu_i} = \phi_i \cdot \frac{1}{\sqrt{2\pi\sigma_i^2}} \cdot \exp\left(-\frac{(x-\mu_i)^2}{2\sigma_i^2}\right) \cdot \frac{x - \mu_i}{\sigma_i^2} \\[5pt]
\frac{\partial G_i}{\partial \mu_i} = G_i(x) \cdot \frac{x - \mu_i}{\sigma_i^2}
$$