Skip to contents

Beyond Cluster Counts: Dual-Anchor Elicitation for Weight Control

JoonHo Lee ()

2026-08-22

Source: vignettes/dual-anchor.Rmd
dual-anchor.Rmd

Overview

The previous vignettes focused on calibrating a Gamma hyperprior for α\alpha based on your expectations about the number of clusters KJK_J. This approach is intuitive: you specify how many groups you expect, and the DPprior package finds the corresponding prior on α\alpha.

However, as Vicentini & Jermyn (2025) have pointed out, matching KJK_J alone can lead to unintended prior behavior on the cluster weights. Specifically, a prior that produces the “right” expected number of clusters may simultaneously imply a large tail probability for the first size-biased cluster mass WSBW_{SB}. This is not the same estimand as the largest cluster mass WmaxW_{max}.

This vignette introduces the dual-anchor framework, which allows you to control both:

  1. Anchor 1: The number of occupied clusters KJK_J (as before)
  2. Anchor 2: A named weight quantity (e.g., WSBW_{SB} or ρ\rho) that captures how mass is distributed across clusters

By the end of this vignette, you will understand:

  • Why K-only calibration can produce surprising weight behavior
  • Two key weight quantities: WSBW_{SB} and ρ\rho (co-clustering probability)
  • How to use DPprior_dual_hard() for a verified inequality
  • How to use DPprior_dual_soft() for an explicit fixed-scale trade-off
  • When and how to choose between hard constraints, soft targets, and different λ\lambda values

1. The Hidden Problem: Unintended Weight Priors

1.1 What K-Only Calibration Misses

Consider a researcher analyzing a multisite educational trial with 50 sites, expecting that treatment effects will cluster into approximately 5 distinct patterns. Using the K-only calibration approach from the Applied Guide, they obtain:

J <- 50
mu_K <- 5

target_K <- DPprior_target_K(J = J, mu_K = mu_K, var_K = 8)
fit_K <- DPprior_fit(
  J = J, target_K = target_K, method = "A2-MN",
  check_diagnostics = FALSE
)
print(fit_K)
#> DPprior Prior Elicitation Result
#> ============================================= 
#> 
#> Schema: dpprior.result/1
#> Method: A2-MN (mode: a2_moment)
#> Status: converged; usable: yes; verified: yes
#> 
#> Target (J = 50):
#>   E[K_J]   = 5.0000
#>   Var(K_J) = 8.0000
#> 
#> Canonical candidate:
#>   alpha ~ Gamma(a = 2.0361, b = 1.6051)
#>   Achieved E[K_J] = 5.000000; Var(K_J) = 8.000000
#>   Maximum absolute moment residual = 1.89e-12
#> 
#> Canonical guidance: scaled component residuals and independent higher-order verification passed

The fit looks reasonable: the prior expects about 5 clusters with appropriate uncertainty. But what does this prior imply about how mass is distributed across those clusters?

# Examine the implied weight behavior
cat("Weight diagnostics for K-only prior:\n")
#> Weight diagnostics for K-only prior:
cat("  E[W_SB] =",
    round(mean_w1(fit_K$parameters$a, fit_K$parameters$b), 3), "\n")
#>   E[W_SB] = 0.501
cat("  P(W_SB > 0.3) =",
    round(prob_wsb_exceeds(0.3, fit_K$parameters$a, fit_K$parameters$b), 3), "\n")
#>   P(W_SB > 0.3) = 0.665
cat("  P(W_SB > 0.5) =",
    round(prob_wsb_exceeds(0.5, fit_K$parameters$a, fit_K$parameters$b), 3), "\n")
#>   P(W_SB > 0.5) = 0.481
cat("  P(W_SB > 0.9) =",
    round(prob_wsb_exceeds(0.9, fit_K$parameters$a, fit_K$parameters$b), 3), "\n")
#>   P(W_SB > 0.9) = 0.163

Here is the key insight: even though we expect only 5 clusters on average, there is a substantial probability (nearly 50%!) that a randomly selected unit belongs to a cluster containing more than half of all units. This is much more concentrated than many researchers would intuitively expect from “5 clusters.”

1.2 Why Does This Happen?

The relationship between KJK_J and the cluster weights is not as tight as one might expect. The expected number of clusters follows approximately:

𝔼[KJ∣α]≈1+αlog⁡J \mathbb{E}[K_J \mid \alpha] \approx 1 + \alpha \log J

while the expected first stick-breaking weight is:

𝔼[w1∣α]=11+α \mathbb{E}[w_1 \mid \alpha] = \frac{1}{1 + \alpha}

For J=50J = 50 and an expected cluster count around 5, we need α≈1\alpha \approx 1, which implies 𝔼[w1∣α]≈0.5\mathbb{E}[w_1 \mid \alpha] \approx 0.5. But because α\alpha is random under our Gamma hyperprior, the tails of the α\alpha distribution can induce extreme weight behavior: small α\alpha values lead to w1w_1 close to 1, while large α\alpha values fragment the weights.

# Illustrate the conditional relationship
alpha_grid <- c(0.5, 1, 2, 5, 10, 20)
w_grid <- seq(0.01, 0.99, length.out = 200)

cond_df <- do.call(rbind, lapply(alpha_grid, function(a) {
  data.frame(
    w = w_grid,
    density = dbeta(w_grid, 1, a),
    alpha = paste0("α = ", a)
  )
}))
cond_df$alpha <- factor(cond_df$alpha, levels = paste0("α = ", alpha_grid))

ggplot(cond_df, aes(x = w, y = density, color = alpha)) +
  geom_line(linewidth = 0.8) +
  scale_color_viridis_d(option = "plasma", end = 0.85) +
  labs(
    x = expression(w[1]),
    y = "Conditional Density",
    title = expression("Conditional Distribution of " * w[1] * " | " * alpha),
    color = NULL
  ) +
  theme_minimal() +
  theme(legend.position = "right")
When α is small, w₁ concentrates near 1; when α is large, w₁ concentrates near 0. A diffuse prior on α spans both extremes.

When α is small, w₁ concentrates near 1; when α is large, w₁ concentrates near 0. A diffuse prior on α spans both extremes.

2. Understanding Weight Distributions

Before diving into dual-anchor calibration, let us examine the two weight anchors available in the DPprior package.

2.1 Stick-Breaking Weights Refresher

Under Sethuraman’s stick-breaking representation of the Dirichlet Process, the random mixing measure GG is:

G=∑h=1∞whδθh G = \sum_{h=1}^{\infty} w_h \, \delta_{\theta_h}

where the weights are constructed via:

vh∼iidBeta(1,α),w1=v1,wh=vh∏ℓ<h(1−vℓ)(h≥2) v_h \overset{iid}{\sim} \text{Beta}(1, \alpha), \quad w_1 = v_1, \quad w_h = v_h \prod_{\ell < h} (1 - v_\ell) \quad (h \geq 2)

The sequence (w1,w2,…)(w_1, w_2, \ldots) follows the GEM(α\alpha) distribution and is in size-biased order, not ranked by magnitude.

2.2 Anchor 2a: The First Stick-Breaking Weight (w1w_1)

The first weight w1w_1 has a particularly tractable distribution. Conditionally, w1∣α∼Beta(1,α)w_1 \mid \alpha \sim \text{Beta}(1, \alpha).

Under the hyperprior α∼Gamma(a,b)\alpha \sim \text{Gamma}(a, b), the marginal distribution of w1w_1 has fully closed-form expressions:

CDF: Fw1(x∣a,b)=1−(bb−log⁡(1−x))a F_{w_1}(x \mid a, b) = 1 - \left( \frac{b}{b - \log(1-x)} \right)^a

Quantile function: Qw1(u∣a,b)=1−exp⁡(b[1−(1−u)−1/a]) Q_{w_1}(u \mid a, b) = 1 - \exp\left( b \left[ 1 - (1-u)^{-1/a} \right] \right)

First-size-biased-weight tail probability: P(w1>t∣a,b)=(bb−log⁡(1−t))a P(w_1 > t \mid a, b) = \left( \frac{b}{b - \log(1-t)} \right)^a

Important interpretability caveat: w1w_1 is in GEM (size-biased) order, not the largest cluster weight. A faithful interpretation is:

“w1w_1 is the asymptotic proportion of the cluster containing a randomly selected unit.”

This is a meaningful named tail diagnostic: if P(WSB>0.5)P(W_{SB}>0.5) is high, then a randomly selected unit is likely to belong to a cluster that contains more than half the population. It is not P(Wmax>0.5)P(W_{max}>0.5).

# Using the closed-form w1 functions
a <- fit_K$parameters$a
b <- fit_K$parameters$b

cat("W_SB distribution under K-only prior:\n")
#> W_SB distribution under K-only prior:
cat("  Mean:     ", round(mean_w1(a, b), 4), "\n")
#>   Mean:      0.5014
cat("  Variance: ", round(var_w1(a, b), 4), "\n")
#>   Variance:  0.1006
cat("  Median:   ", round(quantile_w1(0.5, a, b), 4), "\n")
#>   Median:    0.4784
cat("  90th %ile:", round(quantile_w1(0.9, a, b), 4), "\n")
#>   90th %ile: 0.9655
cat("\nNamed W_SB tail probabilities:\n")
#> 
#> Named W_SB tail probabilities:
cat("  P(W_SB > 0.3):", round(prob_wsb_exceeds(0.3, a, b), 4), "\n")
#>   P(W_SB > 0.3): 0.6646
cat("  P(W_SB > 0.5):", round(prob_wsb_exceeds(0.5, a, b), 4), "\n")
#>   P(W_SB > 0.5): 0.4815
cat("  P(W_SB > 0.9):", round(prob_wsb_exceeds(0.9, a, b), 4), "\n")
#>   P(W_SB > 0.9): 0.1634

2.3 Anchor 2b: Co-Clustering Probability (ρ\rho)

The co-clustering probability is defined as:

ρ=∑h=1∞wh2 \rho = \sum_{h=1}^{\infty} w_h^2

This has a natural interpretation: ρ\rho equals the probability that two randomly selected units belong to the same cluster. This may be more intuitive for applied researchers to elicit.

The conditional moments are simple:

𝔼[ρ∣α]=11+α,Var(ρ∣α)=2α(1+α)2(2+α)(3+α) \mathbb{E}[\rho \mid \alpha] = \frac{1}{1 + \alpha}, \quad \text{Var}(\rho \mid \alpha) = \frac{2\alpha}{(1+\alpha)^2(2+\alpha)(3+\alpha)}

Note that 𝔼[ρ∣α]=𝔼[w1∣α]\mathbb{E}[\rho \mid \alpha] = \mathbb{E}[w_1 \mid \alpha], so their means are equal, but their full distributions differ.

# Using the rho (co-clustering) functions
cat("Co-clustering probability ρ under K-only prior:\n")
#> Co-clustering probability ρ under K-only prior:
cat("  E[ρ]:     ", round(mean_rho(a, b), 4), "\n")
#>   E[ρ]:      0.5014
cat("  Var(ρ):   ", round(var_rho(a, b), 6), "\n")
#>   Var(ρ):    0.064994
cat("  SD(ρ):    ", round(sqrt(var_rho(a, b)), 4), "\n")
#>   SD(ρ):     0.2549

# Compare with w1
cat("\nNote: E[w₁] = E[ρ] =", round(mean_w1(a, b), 4), "\n")
#> 
#> Note: E[w₁] = E[ρ] = 0.5014
cat("But: Var(w₁) =", round(var_w1(a, b), 4), "≠ Var(ρ) =", 
    round(var_rho(a, b), 6), "\n")
#> But: Var(w₁) = 0.1006 ≠ Var(ρ) = 0.064994

3. The Dual-Anchor Framework

3.1 The Core Idea

The v2 dual-anchor framework starts from a usable, verified K-only fit and then selects one of two explicit contracts:

  • Anchor 1: Match your beliefs about KJK_J (number of clusters)
  • Anchor 2: Match a named WSBW_{SB} estimand (tail probability, mean, or quantile), or use the separately certified WmaxW_{max} tail upper bound

Use hard mode when Anchor 2 is a genuine inequality that must be independently verified. It minimizes the fixed-input-scale K loss subject to the named constraint,

mina>0,b>0LK(a,b)subject togw(a,b)≤0. \min_{a>0,b>0} L_K(a,b) \quad\text{subject to}\quad g_w(a,b) \leq 0.

Use soft mode when both anchors are preferences and a trade-off is scientifically appropriate:

(a*,b*)=arg⁡mina>0,b>0[λ⋅LK(a,b)+(1−λ)⋅Lw(a,b)] (a^*, b^*) = \arg\min_{a > 0, b > 0} \left[ \lambda \cdot L_K(a, b) + (1 - \lambda) \cdot L_w(a, b) \right]

where LKL_K measures the discrepancy from the K target and LwL_w measures the discrepancy from the weight target.

Hard mode has no λ\lambda. A hard result is decision-ready only when its mode-specific contract and the result quartet show an approved status with usable = TRUE and verified = TRUE; the independent certificate is in constraint$satisfied.

3.2 The Role of λ\lambda in Soft Mode

The parameter λ∈(0,1]\lambda \in (0, 1] controls the soft trade-off between the two anchors:

λ\lambda Value Interpretation
λ=1\lambda = 1 K-only calibration (ignores weight anchor)
0.5≤λ<10.5 \leq \lambda < 1 K remains primary, weight provides secondary control
0<λ<0.50 < \lambda < 0.5 Weight becomes primary (rare in practice)

In most applications, we recommend λ∈[0.5,0.9]\lambda \in [0.5, 0.9], keeping the cluster count as the primary anchor while using the weight constraint to avoid unintended dominance behavior.

Soft-mode λ\lambda is not a probability and never certifies an inequality. Inspect tradeoff, residuals, and the result quartet before using a soft fit downstream.

3.3 Fixed Input-Derived Scaling

The K and weight losses have different units. DPprior_dual_soft() records fixed scales derived from the input fit and target, then uses those same scales throughout optimization and verification. They are available in fit$tradeoff$scales and fit$computation$scaling. New code does not choose among legacy loss_type heuristics.

4. Using the v2 Dual-Anchor APIs

Begin with DPprior_target_K() and DPprior_fit(). Then choose hard or soft mode from the meaning of the weight statement—not from which optimizer happens to return a convenient answer.

4.1 Basic Usage: Verified Probability Inequality

Suppose you want to ensure that the probability of extreme dominance is limited. Specifically, you want P(w1>0.5)≤0.25P(w_1 > 0.5) \leq 0.25:

# Step 1: construct and fit a canonical K target
target_K <- DPprior_target_K(J = 50, mu_K = 5, var_K = 8)
fit_K <- DPprior_fit(
  J = 50, target_K = target_K, method = "A2-MN",
  check_diagnostics = FALSE
)

cat("K-only prior:\n")
#> K-only prior:
cat("  Gamma(", round(fit_K$parameters$a, 3), ",",
    round(fit_K$parameters$b, 3), ")\n")
#>   Gamma( 2.036 , 1.605 )
cat("  P(w₁ > 0.5) =",
    round(prob_wsb_exceeds(
      0.5, fit_K$parameters$a, fit_K$parameters$b
    ), 3), "\n\n")
#>   P(w₁ > 0.5) = 0.481

# Step 2: impose and independently verify an at-most constraint
fit_dual <- DPprior_dual_hard(
  fit = fit_K,
  constraint = list(
    metric = "wsb_tail", threshold = 0.5,
    relation = "<=", bound = 0.25
  )
)

cat("Hard dual-anchor prior:\n")
#> Hard dual-anchor prior:
cat("  Gamma(", round(fit_dual$parameters$a, 3), ",",
    round(fit_dual$parameters$b, 3), ")\n")
#>   Gamma( 8.004 , 3.666 )
cat("  P(w₁ > 0.5) =", round(fit_dual$achieved$weight$value, 3), "\n")
#>   P(w₁ > 0.5) = 0.25
cat("  Constraint verified:", fit_dual$constraint$satisfied, "\n")
#>   Constraint verified: TRUE
fit_dual[c("status", "usable", "verified", "message")]
#> $status
#> [1] "boundary"
#> 
#> $usable
#> [1] TRUE
#> 
#> $verified
#> [1] TRUE
#> 
#> $message
#> [1] "Hard inequality solution passed verification with the constraint active."

Let us visualize the effect on both the KJK_J and w1w_1 distributions:

# Prepare data for visualization
logS <- compute_log_stirling(J)

# K distributions
k_grid <- 1:20
k_df <- data.frame(
  K = rep(k_grid, 2),
  probability = c(
    pmf_K_marginal(
      J, fit_K$parameters$a, fit_K$parameters$b, logS = logS
    )[k_grid],
    pmf_K_marginal(
      J, fit_dual$parameters$a, fit_dual$parameters$b, logS = logS
    )[k_grid]
  ),
  Prior = rep(c("K-only", "Dual-anchor"), each = length(k_grid))
)

# w1 distributions
w_grid <- seq(0.01, 0.95, length.out = 200)
w_df <- data.frame(
  w1 = rep(w_grid, 2),
  density = c(
    density_w1(w_grid, fit_K$parameters$a, fit_K$parameters$b),
    density_w1(w_grid, fit_dual$parameters$a, fit_dual$parameters$b)
  ),
  Prior = rep(c("K-only", "Dual-anchor"), each = length(w_grid))
)

# Create plots
p1 <- ggplot(k_df, aes(x = K, y = probability, color = Prior)) +
  geom_point(size = 1.5) +
  geom_line(linewidth = 0.8) +
  scale_color_manual(values = palette_2) +
  labs(x = "Number of Clusters (K)", y = "Probability",
       title = "Marginal PMF of K") +
  theme_minimal() +
  theme(legend.position = "bottom")

p2 <- ggplot(w_df, aes(x = w1, y = density, color = Prior)) +
  geom_line(linewidth = 1) +
  geom_vline(xintercept = 0.5, linetype = "dashed", color = "gray50") +
  scale_color_manual(values = palette_2) +
  labs(x = expression(w[1]), y = "Density",
       title = expression("Marginal Density of " * w[1])) +
  theme_minimal() +
  theme(legend.position = "bottom")

gridExtra::grid.arrange(p1, p2, ncol = 2)
Comparison of K-only vs dual-anchor priors: the dual-anchor approach shifts mass away from extreme w₁ values while maintaining similar K expectations.

Comparison of K-only vs dual-anchor priors: the dual-anchor approach shifts mass away from extreme w₁ values while maintaining similar K expectations.

4.2 Comparing Hard and Soft Semantics

The same numeric value can express different scientific claims. Here the hard fit certifies an upper bound, while the soft fit treats 0.25 as an equality preference with λ=0.5\lambda=0.5:

soft_target <- list(
  metric = "wsb_tail", threshold = 0.5,
  relation = "target", value = 0.25
)
fit_soft <- DPprior_dual_soft(
  fit = fit_K, target = soft_target, lambda = 0.5
)

fits_compare <- list(K_only = fit_K, Hard = fit_dual, Soft = fit_soft)
hard_soft_compare <- data.frame(
  Method = names(fits_compare),
  a = vapply(fits_compare, function(x) x$parameters$a, numeric(1)),
  b = vapply(fits_compare, function(x) x$parameters$b, numeric(1)),
  E_K = vapply(fits_compare, function(x) x$achieved$K$mean, numeric(1)),
  P_w1_gt_50 = vapply(fits_compare, function(x) {
    prob_wsb_exceeds(0.5, x$parameters$a, x$parameters$b)
  }, numeric(1)),
  Status = vapply(fits_compare, function(x) x$status, character(1)),
  Usable = vapply(fits_compare, function(x) x$usable, logical(1)),
  Verified = vapply(fits_compare, function(x) x$verified, logical(1))
)

knitr::kable(
  hard_soft_compare,
  digits = 3,
  caption = paste(
    "Hard inequality versus fixed-scale soft target for",
    "P(w₁ > 0.5) = 0.25"
  )
)
Hard inequality versus fixed-scale soft target for P(w₁ > 0.5) = 0.25
Method a b E_K P_w1_gt_50 Status Usable Verified
K_only K_only 2.036 1.605 5.000 0.481 converged TRUE TRUE
Hard Hard 8.004 3.666 7.329 0.250 boundary TRUE TRUE
Soft Soft 2.735 1.898 5.480 0.427 converged TRUE TRUE

Only the hard result has constraint$satisfied; a soft result instead records component and total losses in tradeoff.

4.3 Alternative Hard Constraint: Quantile

You can also specify a quantile constraint. For example, to ensure that the 90th percentile of w1w_1 does not exceed 0.6:

fit_dual_q <- DPprior_dual_hard(
  fit = fit_K,
  constraint = list(
    metric = "wsb_quantile", probability = 0.9,
    relation = "<=", bound = 0.6
  )
)

cat("Hard dual-anchor (quantile constraint):\n")
#> Hard dual-anchor (quantile constraint):
cat("  Gamma(", round(fit_dual_q$parameters$a, 3), ",",
    round(fit_dual_q$parameters$b, 3), ")\n")
#>   Gamma( 21.319 , 8.034 )
cat("  90th percentile of w₁:",
    round(fit_dual_q$achieved$weight$value, 3), "\n")
#>   90th percentile of w₁: 0.6
cat("  Constraint verified:", fit_dual_q$constraint$satisfied, "\n")
#>   Constraint verified: TRUE

4.4 Alternative Soft Target: Mean

For constraints based on 𝔼[w1]\mathbb{E}[w_1]:

fit_dual_m <- DPprior_dual_soft(
  fit = fit_K,
  target = list(
    metric = "wsb_mean", relation = "target", value = 0.35
  ),
  lambda = 0.6
)

cat("Soft dual-anchor (mean target):\n")
#> Soft dual-anchor (mean target):
cat("  Gamma(", round(fit_dual_m$parameters$a, 3), ",",
    round(fit_dual_m$parameters$b, 3), ")\n")
#>   Gamma( 2.309 , 1.726 )
cat("  E[w₁]:", round(fit_dual_m$achieved$weight$value, 3), "\n")
#>   E[w₁]: 0.483
cat("  (Target was 0.35)\n")
#>   (Target was 0.35)
print(fit_dual_m$tradeoff[c("lambda", "K_loss", "weight_loss", "total_loss")])
#> $lambda
#> [1] 0.6
#> 
#> $K_loss
#> [1] 0.001590141
#> 
#> $weight_loss
#> [1] 0.01770233
#> 
#> $total_loss
#> [1] 0.008035014
fit_dual_m[c("status", "usable", "verified", "message")]
#> $status
#> [1] "converged"
#> 
#> $usable
#> [1] TRUE
#> 
#> $verified
#> [1] TRUE
#> 
#> $message
#> [1] "soft trade-off optimizer and independent verification passed"

5. Elicitation Questions for Weight Anchors

Translating substantive knowledge into weight constraints requires careful framing. Here are suggested elicitation questions:

For w1w_1 (Size-Biased Weight)

These questions focus on what happens when you randomly sample a unit:

“If you randomly select a site from your study, what proportion of all sites would you expect to share its effect pattern? What’s a reasonable median estimate?”

“Is there a significant chance (say, > 30%) that a randomly selected site belongs to an effect group that contains more than half of all sites?”

For ρ\rho (Co-Clustering Probability)

These questions may be more intuitive for many researchers:

“If you pick two sites at random, what’s the probability that they belong to the same effect group?”

“How likely is it that any two randomly chosen sites would show substantively similar treatment effects?”

Translation Examples

Elicited Belief Formal Constraint
“Randomly selected site’s group is at most 25% (median)” Hard: metric="wsb_quantile", probability=.5, relation="<=", bound=.25
“At most 20% chance that a random site is in a dominant (>50%) group” Hard: metric="wsb_tail", threshold=.5, relation="<=", bound=.2
“Two random sites have about 20% chance of the same group” Diagnose ρ\rho separately; do not relabel it as a v2 WSBW_{SB} target

6. Exploring the Trade-off: λ\lambda Sensitivity

Different values of λ\lambda produce different trade-offs between matching the K anchor and the weight anchor. Let us explore this systematically:

lambda_grid <- c(1.0, 0.9, 0.7, 0.5, 0.3)
w1_target <- list(
  metric = "wsb_tail", threshold = 0.5,
  relation = "target", value = 0.2
)

lambda_results <- lapply(lambda_grid, function(lam) {
  fit <- DPprior_dual_soft(
    fit = fit_K, target = w1_target, lambda = lam
  )
  
  list(
    lambda = lam,
    a = fit$parameters$a,
    b = fit$parameters$b,
    E_K = fit$achieved$K$mean,
    P_w1_gt_50 = prob_wsb_exceeds(
      0.5, fit$parameters$a, fit$parameters$b
    ),
    E_w1 = mean_w1(fit$parameters$a, fit$parameters$b),
    status = fit$status,
    usable = fit$usable,
    verified = fit$verified
  )
})

# Create comparison table
lambda_df <- do.call(rbind, lapply(lambda_results, function(r) {
  data.frame(
    lambda = r$lambda,
    a = round(r$a, 3),
    b = round(r$b, 3),
    E_K = round(r$E_K, 2),
    P_w1_gt_50 = round(r$P_w1_gt_50, 3),
    E_w1 = round(r$E_w1, 3),
    status = r$status,
    usable = r$usable,
    verified = r$verified
  )
}))

knitr::kable(
  lambda_df,
  col.names = c(
    "λ", "a", "b", "E[K]", "P(w₁>0.5)", "E[w₁]",
    "Status", "Usable", "Verified"
  ),
  caption = paste(
    "Effect of λ under the fixed-scale soft contract",
    "(target: P(w₁>0.5) = 0.2)"
  )
)
Effect of λ under the fixed-scale soft contract (target: P(w₁>0.5) = 0.2)
λ a b E[K] P(w₁>0.5) E[w₁] Status Usable Verified
1.0 2.036 1.605 5.00 0.481 0.501 converged TRUE TRUE
0.9 2.150 1.655 5.09 0.471 0.493 converged TRUE TRUE
0.7 2.447 1.781 5.30 0.447 0.474 converged TRUE TRUE
0.5 2.899 1.962 5.58 0.416 0.450 converged TRUE TRUE
0.3 3.708 2.263 6.00 0.371 0.417 converged TRUE TRUE 
# Visualize the trade-off
tradeoff_df <- data.frame(
  lambda = sapply(lambda_results, `[[`, "lambda"),
  E_K = sapply(lambda_results, `[[`, "E_K"),
  P_w1 = sapply(lambda_results, `[[`, "P_w1_gt_50")
)

ggplot(tradeoff_df, aes(x = E_K, y = P_w1)) +
  geom_path(color = "gray50", linewidth = 0.8) +
  geom_point(aes(color = factor(lambda)), size = 4) +
  geom_hline(yintercept = 0.2, linetype = "dashed", color = "#E41A1C", alpha = 0.7) +
  geom_vline(xintercept = 5, linetype = "dashed", color = "#377EB8", alpha = 0.7) +
  annotate("text", x = 4.7, y = 0.22, label = "w₁ target", hjust = 1, color = "#E41A1C") +
  annotate("text", x = 5.1, y = 0.45, label = "K target", hjust = 0, color = "#377EB8") +
  scale_color_viridis_d(option = "viridis", begin = 0.2, end = 0.8) +
  labs(
    x = "E[K] (Cluster Count)",
    y = expression(P(w[1] > 0.5)),
    title = "Trade-off Between K and Weight Anchors",
    color = "λ"
  ) +
  theme_minimal() +
  theme(legend.position = "right")
Trade-off curve: reducing λ brings P(w₁ > 0.5) closer to target at the cost of K deviation.

Trade-off curve: reducing λ brings P(w₁ > 0.5) closer to target at the cost of K deviation.

7. Practical Recommendations

7.1 When to Use Dual-Anchor Calibration

Consider dual-anchor calibration when:

  • A predeclared WSBW_{SB} policy is triggered: For example, the reviewed fit has P(WSB>0.5)P(W_{SB}>0.5) above an application-specific action threshold

  • Strong prior belief against concentration: You believe effects should be reasonably spread across groups, not concentrated

  • Low-information settings: When the data may not strongly update the prior, unintended weight behavior can persist to the posterior

Do not use dual-anchor when:

  • You are comfortable with the weight implications of your K-only prior
  • You have strong data that will overwhelm any prior specification
  • The additional complexity is not justified for your application

7.2 Choosing Hard or Soft Mode

Scientific statement API Evidence to inspect
A weight inequality must hold DPprior_dual_hard() constraint$satisfied, residual, tolerance, verification
K and weight targets are competing preferences DPprior_dual_soft() tradeoff, residuals, result quartet

Do not report a soft target as if it were a hard constraint. Conversely, do not introduce λ\lambda into a hard inequality.

7.3 Choosing λ\lambda in Soft Mode

Situation Recommended λ\lambda
K anchor is well-justified, weight is secondary constraint 0.8 – 0.9
Both anchors are equally important 0.5 – 0.7
Weight behavior is primary concern 0.3 – 0.5

In practice, we recommend starting with λ=0.7\lambda = 0.7 and adjusting based on how closely each target is achieved.

7.4 Choosing a Supported Weight Estimand

Estimand Advantages Best For
WSBW_{SB} tail Closed form; explicitly size-biased Probability inequalities or targets
WSBW_{SB} mean Simple scalar interpretation Mean preference or bound
WSBW_{SB} quantile Closed form Quantile preference or bound
certified WmaxW_{max} tail upper bound Conservative certified evidence At-most safety constraints

The co-clustering quantity ρ\rho remains a useful diagnostic, and 𝔼[w1]=𝔼[ρ]\mathbb{E}[w_1]=\mathbb{E}[\rho], but the v2 hard/soft public target grammar does not silently relabel a WSBW_{SB} claim as a direct ρ\rho claim.

8. When Dual-Anchor Constraints Are Infeasible

Not all combinations of K and weight constraints are achievable. Hard mode does not return a failed inequality as if it were satisfied. It raises a typed condition whose result retains the complete canonical evidence:

# An aggressive weight target that conflicts with K target
hard_condition <- NULL
fit_aggressive <- tryCatch(
  DPprior_dual_hard(
    fit = fit_K,
    constraint = list(
      metric = "wsb_tail", threshold = 0.5,
      relation = "<=", bound = 0.05
    )
  ),
  error = function(condition) {
    hard_condition <<- condition
    condition$result
  }
)

cat("Aggressive dual-anchor attempt:\n")
#> Aggressive dual-anchor attempt:
cat("  Target: P(w₁ > 0.5) = 0.05\n")
#>   Target: P(w₁ > 0.5) = 0.05
fit_aggressive[c("status", "usable", "verified", "message")]
#> $status
#> [1] "boundary"
#> 
#> $usable
#> [1] TRUE
#> 
#> $verified
#> [1] TRUE
#> 
#> $message
#> [1] "Hard inequality solution passed verification at a declared log-parameter boundary."
if (!is.null(fit_aggressive$parameters)) {
  cat("  Achieved: P(w₁ > 0.5) =",
      round(fit_aggressive$achieved$weight$value, 3), "\n")
  cat("  Constraint satisfied:",
      fit_aggressive$constraint$satisfied, "\n")
}
#>   Achieved: P(w₁ > 0.5) = 0.05 
#>   Constraint satisfied: TRUE
if (!is.null(hard_condition)) {
  cat("  Typed condition:", class(hard_condition)[1], "\n")
}

When constraints are too strict, consider:

  1. Relaxing one of the targets: Accept a higher named WSBW_{SB} tail probability or a different K expectation

  2. Changing the scientific contract: If compromise is defensible, use soft mode and perform an explicit λ\lambda sensitivity analysis

  3. Using a single anchor: Sometimes K-only or weight-only calibration is more appropriate

Summary

Concept Key Point
K-only calibration May induce unintended weight behavior
w1w_1 (first weight) Size-biased cluster mass; closed-form distribution
ρ\rho (co-clustering) P(two random units share a cluster); intuitive
Hard dual anchor Verified named inequality; no λ\lambda
Soft dual anchor Fixed-scale trade-off between K and weight fit via λ\lambda
Recommended workflow 1. DPprior_target_K() → 2. DPprior_fit() → 3. hard or soft dual API

Legacy v2.x Compatibility

DPprior_dual() is retained throughout v2.x for existing equality-loss workflows. It is an approximate legacy adapter, not a hard-constraint certificate: its canonical result is visibly status = "approximate", usable = TRUE, and verified = FALSE. It will not be removed before v3.0 and a migration review. New analyses should use DPprior_dual_hard() or DPprior_dual_soft().

# Compatibility only; do not interpret this as a verified inequality.
legacy_fit <- DPprior_dual(
  fit = fit_K,
  w1_target = list(prob = list(threshold = 0.5, value = 0.25)),
  lambda = 0.7,
  loss_type = "adaptive"
)
legacy_fit[c("status", "usable", "verified", "message")]

What’s Next?

References

Vicentini, S., & Jermyn, I. H. (2025). Prior selection for the precision parameter of Dirichlet process mixture models. arXiv:2502.00864. https://doi.org/10.48550/arXiv.2502.00864

Sethuraman, J. (1994). A constructive definition of Dirichlet priors. Statistica Sinica, 4(2), 639–650.

For questions or feedback, please visit the GitHub repository.