Skip to content

API Reference

stdpy — Python measures

The stdpy sub-module is a full NumPy/SciPy implementation of every measure.

Entry point

compute_all_measures(stim, resp, conf, n_ratings, *, matlab_compat=False, return_type='array')

Compute all 26 meta-signal measures mimicking the MATLAB output.

Returns 26 values corresponding to (see MEASURE_NAMES for the exact, ordered list): [ 1: meta_d, 2: AUC2, 3: gamma, 4: phi, 5: deltaConf, 6: M_ratio, 7: AUC2_ratio, 8: gamma_ratio, 9: phi_ratio, 10: deltaConf_ratio, 11: M_diff, 12: AUC2_diff, 13: gamma_diff, 14: phi_diff, 15: deltaConf_diff, 16: metaNoise, 17: metaUncertainty, 18: dprime, 19: c, 20: mean_conf, 21: logL, 22: AIC, 23: BIC, 24: AICc, 25: k, 26: n ]

matlab_compat

Forwarded to fit_meta_d_mle / compute_meta_uncertainty for single-start / high-budget MATLAB-like optimizer behaviour.

return_type : 'array' (default) — numpy.ndarray of 26 values, positional (backward compatible). 'dict' — {measure_name: value}. 'dataframe' — single-row pandas.DataFrame with measure names as columns.

Source code in src/metasignal/stdpy/compute_all.py 
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
def compute_all_measures(
    stim: np.ndarray,
    resp: np.ndarray,
    conf: np.ndarray,
    n_ratings: int,
    *,
    matlab_compat: bool = False,
    return_type: str = "array",
) -> "np.ndarray | dict | pd.DataFrame":
    """Compute all 26 meta-signal measures mimicking the MATLAB output.

    Returns 26 values corresponding to (see ``MEASURE_NAMES`` for the exact,
    ordered list):
    [
      1: meta_d, 2: AUC2, 3: gamma, 4: phi, 5: deltaConf,
      6: M_ratio, 7: AUC2_ratio, 8: gamma_ratio, 9: phi_ratio, 10: deltaConf_ratio,
      11: M_diff, 12: AUC2_diff, 13: gamma_diff, 14: phi_diff, 15: deltaConf_diff,
      16: metaNoise, 17: metaUncertainty, 18: dprime, 19: c, 20: mean_conf,
      21: logL, 22: AIC, 23: BIC, 24: AICc, 25: k, 26: n
    ]

    matlab_compat :
        Forwarded to ``fit_meta_d_mle`` / ``compute_meta_uncertainty`` for
        single-start / high-budget MATLAB-like optimizer behaviour.
    return_type :
        ``'array'`` *(default)* — ``numpy.ndarray`` of 26 values, positional
        (backward compatible).
        ``'dict'`` — ``{measure_name: value}``.
        ``'dataframe'`` — single-row ``pandas.DataFrame`` with measure names
        as columns.
    """
    if return_type not in _RETURN_TYPES:
        raise ValueError(
            f"return_type must be one of {sorted(_RETURN_TYPES)}, got {return_type!r}"
        )
    stim = np.asarray(stim, dtype=float)
    resp = np.asarray(resp, dtype=float)
    conf = np.asarray(conf, dtype=float)

    # Remove NaNs
    valid = ~np.isnan(stim) & ~np.isnan(resp) & ~np.isnan(conf)
    stim = stim[valid]
    resp = resp[valid]
    conf = conf[valid]

    if len(stim) == 0:
        return _format_result(np.full(26, np.nan), return_type)

    if len(np.unique(stim)) > 2:
        msg = f"stim must be binary (2 classes); found {len(np.unique(stim))} distinct values"
        raise ValueError(msg)
    if len(np.unique(resp)) > 2:
        msg = f"resp must be binary (2 classes); found {len(np.unique(resp))} distinct values"
        raise ValueError(msg)

    # Make input 0/1
    stim_max = np.max(stim)
    resp_max = np.max(resp)

    stim_bin = np.zeros_like(stim, dtype=int)
    stim_bin[stim == stim_max] = 1

    resp_bin = np.zeros_like(resp, dtype=int)
    resp_bin[resp == resp_max] = 1

    # Basic quantities
    dprime, c, _ = compute_sdt_resp(stim_bin, resp_bin)
    mean_conf = np.mean(conf)

    if np.array_equal(stim_bin, resp_bin) or dprime == 0 or len(np.unique(conf)) == 1:
        return _format_result(np.full(26, np.nan), return_type)

    # Convert to counts (raw, unpadded — used for all descriptive Type-2 measures)
    nr_s1, nr_s2 = trials_to_counts(stim_bin, resp_bin, conf.astype(int), n_ratings)
    nr_s1 = np.array(nr_s1)
    nr_s2 = np.array(nr_s2)

    # Zero-cell padding (matches MATLAB type2_SDT_MLE.m behaviour):
    # add 1/(2*nRatings) to ALL cells when ANY cell is zero. This padding
    # only stabilizes the meta-d' MLE log-likelihood (which takes logs of
    # cell probabilities); AUC2/Gamma/Phi/DeltaConf below are purely
    # descriptive statistics of the raw counts and must not be smoothed.
    nr_s1_mle = nr_s1
    nr_s2_mle = nr_s2
    if np.any(nr_s1 == 0) or np.any(nr_s2 == 0):
        pad = 1.0 / (2 * n_ratings)
        nr_s1_mle = nr_s1 + pad
        nr_s2_mle = nr_s2 + pad

    # meta-d', M-Ratio, M-Diff, model fit stats
    try:
        meta_d_res = fit_meta_d_mle(
            nr_s1_mle, nr_s2_mle, matlab_compat=matlab_compat
        )
        if not meta_d_res["success"]:
            warnings.warn(
                "fit_meta_d_mle did not converge; meta-d’-derived measures set to NaN.",
                stacklevel=2,
            )
            raise RuntimeError("fit_meta_d_mle failed to converge")
        meta_d = meta_d_res["meta_da"]
        m_ratio = meta_d_res["M_ratio"]
        m_diff  = meta_d_res["M_diff"]
        logL    = meta_d_res["logL"]
        aic     = meta_d_res["AIC"]
        bic     = meta_d_res["BIC"]
        aicc    = meta_d_res["AICc"]
        k_fit   = meta_d_res["k"]
        n_fit   = meta_d_res["n"]
    except (ValueError, RuntimeError):
        meta_d, m_ratio, m_diff = np.nan, np.nan, np.nan
        logL, aic, bic, aicc, k_fit, n_fit = np.nan, np.nan, np.nan, np.nan, np.nan, np.nan

    if dprime < 0.2:  # avoid dividing by a negative or very small d' (matches MATLAB)
        m_ratio = np.nan

    # SDT Expectations for ratios/diffs
    sdt_exp = sdt_expect_conf(nr_s1, nr_s2)
    nr_s1_exp = np.array(sdt_exp["nR_S1_exp"])
    nr_s2_exp = np.array(sdt_exp["nR_S2_exp"])

    # AUC2
    auc2 = compute_type2_auc(nr_s1, nr_s2)
    auc2_exp = compute_type2_auc(nr_s1_exp, nr_s2_exp)
    auc2_ratio = auc2 / auc2_exp if auc2_exp != 0 else np.nan
    auc2_diff = auc2 - auc2_exp

    # Gamma
    gamma = compute_gamma(nr_s1, nr_s2)
    gamma_exp = compute_gamma(nr_s1_exp, nr_s2_exp)
    gamma_ratio = gamma / gamma_exp if gamma_exp != 0 else np.nan
    gamma_diff = gamma - gamma_exp

    # Phi
    phi = compute_phi(nr_s1, nr_s2)
    phi_exp = compute_phi(nr_s1_exp, nr_s2_exp)
    phi_ratio = phi / phi_exp if phi_exp != 0 else np.nan
    phi_diff = phi - phi_exp

    # DeltaConf
    dc_res = compute_delta_conf(nr_s1, nr_s2)
    delta_conf = dc_res["delta_conf"]
    delta_conf_ratio = dc_res["delta_conf_ratio"]
    delta_conf_diff = dc_res["delta_conf_diff"]

    # metaNoise
    try:
        mn_res = compute_meta_noise(stim_bin, resp_bin, conf.astype(int), n_ratings)
        meta_noise = mn_res["meta_noise"]
    except (ValueError, RuntimeError):
        meta_noise = np.nan

    # metaUncertainty
    try:
        meta_uncert = compute_meta_uncertainty(
            stim_bin,
            resp_bin,
            conf.astype(int),
            n_ratings,
            matlab_compat=matlab_compat,
        )
    except (ValueError, RuntimeError):
        meta_uncert = np.nan

    # Combine all
    values = np.array([
        meta_d, auc2, gamma, phi, delta_conf,
        m_ratio, auc2_ratio, gamma_ratio, phi_ratio, delta_conf_ratio,
        m_diff, auc2_diff, gamma_diff, phi_diff, delta_conf_diff,
        meta_noise, meta_uncert, float(dprime), float(c), float(mean_conf),
        logL, aic, bic, aicc, float(k_fit), float(n_fit),
    ])
    return _format_result(values, return_type)

Signal Detection Theory

compute_sdt_resp(stimulus, response)

Compute SDT parameters (d', c, ln_beta).

Parameters:

Name Type Description Default
stimulus ndarray

Vector of 2 values (lower is noise, higher is stimulus).

 required
response ndarray

Vector with the same values as stimulus.

 required

Returns:

Name Type Description
tuple tuple[float, float, float]

(dprime, c, ln_beta)
Source code in src/metasignal/stdpy/core.py 
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
def compute_sdt_resp(
    stimulus: np.ndarray, response: np.ndarray
) -> tuple[float, float, float]:
    """Compute SDT parameters (d', c, ln_beta).

    Args:
        stimulus: Vector of 2 values (lower is noise, higher is stimulus).
        response: Vector with the same values as stimulus.

    Returns:
        tuple: (dprime, c, ln_beta)
    """
    s_min = np.min(stimulus)
    s_max = np.max(stimulus)

    n_s1 = np.sum(stimulus == s_min)
    n_s2 = np.sum(stimulus == s_max)

    if n_s1 == 0 or n_s2 == 0:
        raise ValueError(
            "stimulus must contain both classes; found only one class in input."
        )

    # Determine hit and FA rate
    hit_rate = np.sum((stimulus == s_max) & (response == s_max)) / n_s2
    fa_rate = np.sum((stimulus == s_min) & (response == s_max)) / n_s1

    # Correct for values of 0 or 1
    if hit_rate == 0:
        hit_rate = 0.5 / n_s2
    elif hit_rate == 1:
        hit_rate = 1 - 0.5 / n_s2

    if fa_rate == 0:
        fa_rate = 0.5 / n_s1
    elif fa_rate == 1:
        fa_rate = 1 - 0.5 / n_s1

    # Compute d' and criterion c
    zh = norm.ppf(hit_rate)
    zfa = norm.ppf(fa_rate)

    dprime = zh - zfa
    c = -0.5 * (zh + zfa)
    ln_beta = dprime * c

    return dprime, c, ln_beta

trials_to_counts(stim_id, response, rating, n_ratings)

Convert trial data to response counts for each stimulus category.

Parameters:

Name Type Description Default
stim_id ndarray

0=S1, 1=S2.

 required
response ndarray

0=S1, 1=S2.

 required
rating ndarray

1 to n_ratings.

 required
n_ratings int

Total number of ratings.

 required

Returns:

Name Type Description
tuple tuple[ndarray, ndarray]

(nr_s1, nr_s2)
Source code in src/metasignal/stdpy/core.py 
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
def trials_to_counts(
    stim_id: np.ndarray, response: np.ndarray, rating: np.ndarray, n_ratings: int
) -> tuple[np.ndarray, np.ndarray]:
    """Convert trial data to response counts for each stimulus category.

    Args:
        stim_id: 0=S1, 1=S2.
        response: 0=S1, 1=S2.
        rating: 1 to n_ratings.
        n_ratings: Total number of ratings.

    Returns:
        tuple: (nr_s1, nr_s2)
    """
    # Filter bad trials
    f = (
        (np.isin(stim_id, [0, 1]))
        & (np.isin(response, [0, 1]))
        & (rating >= 1)
        & (rating <= n_ratings)
    )
    n_dropped = int(np.sum(~f))
    if n_dropped > 0:
        warnings.warn(
            f"trials_to_counts: dropped {n_dropped} of {len(stim_id)} trials with "
            f"stim/response outside {{0, 1}} or rating outside [1, {n_ratings}] "
            "(check that n_ratings matches the actual confidence scale).",
            stacklevel=2,
        )
    stim_id_f = stim_id[f]
    response_f = response[f]
    rating_f = rating[f]

    nr_s1 = np.zeros(2 * n_ratings)
    nr_s2 = np.zeros(2 * n_ratings)

    # get tallies of "S1" rating responses for S1 and S2 stim
    for i in range(1, n_ratings + 1):
        nr_s1[i - 1] = np.sum(
            (stim_id_f == 0) & (response_f == 0) & (rating_f == n_ratings + 1 - i)
        )
        nr_s2[i - 1] = np.sum(
            (stim_id_f == 1) & (response_f == 0) & (rating_f == n_ratings + 1 - i)
        )

    # get tallies of "S2" rating responses for S1 and S2 stim
    for i in range(1, n_ratings + 1):
        nr_s1[i + n_ratings - 1] = np.sum(
            (stim_id_f == 0) & (response_f == 1) & (rating_f == i)
        )
        nr_s2[i + n_ratings - 1] = np.sum(
            (stim_id_f == 1) & (response_f == 1) & (rating_f == i)
        )

    return nr_s1, nr_s2

Meta-d' (MLE)

fit_meta_d_mle(nr_s1, nr_s2, s=1.0, *, matlab_compat=False)

Fit meta-d' using Maximum Likelihood Estimation.

Parameters:

Name Type Description Default
nr_s1 ndarray

Counts for S1 stimulus.

 required
nr_s2 ndarray

Counts for S2 stimulus.

 required
s float

Ratio of standard deviations (sd(S1)/sd(S2)). Default is 1.0.

 1.0
matlab_compat bool

If True, use MATLAB fit_meta_d_MLE-like optimizer budgets (maxiter/maxfun = 1e5) and still return the fit even when SciPy reports non-convergence (MATLAB packages the returned fmincon point regardless).

 False

Returns:

Name Type Description
dict dict[str, Any]

Results of the fit.
Source code in src/metasignal/stdpy/metad.py 
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
def fit_meta_d_mle(
    nr_s1: np.ndarray,
    nr_s2: np.ndarray,
    s: float = 1.0,
    *,
    matlab_compat: bool = False,
) -> dict[str, Any]:
    """Fit meta-d' using Maximum Likelihood Estimation.

    Args:
        nr_s1: Counts for S1 stimulus.
        nr_s2: Counts for S2 stimulus.
        s: Ratio of standard deviations (sd(S1)/sd(S2)). Default is 1.0.
        matlab_compat: If True, use MATLAB ``fit_meta_d_MLE``-like optimizer
            budgets (``maxiter``/``maxfun`` = 1e5) and still return the fit
            even when SciPy reports non-convergence (MATLAB packages the
            returned ``fmincon`` point regardless).

    Returns:
        dict: Results of the fit.
    """
    if len(nr_s1) % 2 != 0:
        val_err = "input arrays must have an even number of elements"
        raise ValueError(val_err)
    if len(nr_s1) != len(nr_s2):
        val_err = "input arrays must have the same number of elements"
        raise ValueError(val_err)
    if np.sum(nr_s1) == 0 or np.sum(nr_s2) == 0:
        val_err = "nr_s1 and nr_s2 must each contain at least one trial"
        raise ValueError(val_err)

    n_ratings = len(nr_s1) // 2
    n_criteria = 2 * n_ratings - 1

    # Calculate initial d1 and t1c
    rating_hr = np.cumsum(nr_s2[::-1]) / np.sum(nr_s2)
    rating_far = np.cumsum(nr_s1[::-1]) / np.sum(nr_s1)
    rating_hr = rating_hr[:-1][::-1]
    rating_far = rating_far[:-1][::-1]

    t1_index = n_ratings - 1

    eps = 1e-5
    rating_hr = np.clip(rating_hr, eps, 1 - eps)
    rating_far = np.clip(rating_far, eps, 1 - eps)

    d1 = (1 / s) * norm.ppf(rating_hr[t1_index]) - norm.ppf(rating_far[t1_index])
    t1c1 = (-1 / (1 + s)) * (
        norm.ppf(rating_hr[t1_index]) + norm.ppf(rating_far[t1_index])
    )

    constant_criterion_val = d1 * (t1c1 / d1) if d1 != 0 else 0

    c1 = (-1 / (1 + s)) * (norm.ppf(rating_hr) + norm.ppf(rating_far))
    t2c1_guess = np.delete(c1, t1_index)

    guess = np.concatenate([[d1], t2c1_guess - constant_criterion_val])

    lb = np.concatenate(
        [[-10], -20 * np.ones((n_criteria - 1) // 2), np.zeros((n_criteria - 1) // 2)]
    )
    ub = np.concatenate(
        [[10], np.zeros((n_criteria - 1) // 2), 20 * np.ones((n_criteria - 1) // 2)]
    )
    bounds = list(zip(lb, ub))

    def monotonicity_constraint(params: np.ndarray) -> np.ndarray:
        t2c_rel = params[1:]
        half = (n_criteria - 1) // 2
        full_c = np.concatenate([t2c_rel[:half], [0], t2c_rel[half:]])
        return np.diff(full_c) - 1e-5

    cons = {"type": "ineq", "fun": monotonicity_constraint}

    def neg_log_likelihood(params: np.ndarray) -> float:
        meta_d1_fit = params[0]
        t2c_rel = params[1:]

        half = (n_criteria - 1) // 2
        curr_constant_criterion = meta_d1_fit * (t1c1 / d1) if d1 != 0 else 0

        s1mu = -meta_d1_fit / 2 - curr_constant_criterion
        s1sd = 1.0
        s2mu = meta_d1_fit / 2 - curr_constant_criterion
        s2sd = s1sd / s

        t2c_full = np.concatenate([t2c_rel[:half], [0], t2c_rel[half:]])
        t2c_full_ext = np.concatenate([[-np.inf], t2c_full, [np.inf]])

        nc_rs1 = nr_s1[:n_ratings]
        ni_rs1 = nr_s2[:n_ratings]
        nc_rs2 = nr_s2[n_ratings:]
        ni_rs2 = nr_s1[n_ratings:]

        c_area_rs1 = norm.cdf(0, s1mu, s1sd)
        i_area_rs1 = norm.cdf(0, s2mu, s2sd)
        c_area_rs2 = 1 - norm.cdf(0, s2mu, s2sd)
        i_area_rs2 = 1 - norm.cdf(0, s1mu, s1sd)

        c_area_rs1 = max(c_area_rs1, 1e-10)
        i_area_rs1 = max(i_area_rs1, 1e-10)
        c_area_rs2 = max(c_area_rs2, 1e-10)
        i_area_rs2 = max(i_area_rs2, 1e-10)

        log_l = 0.0
        for i in range(n_ratings):
            pc = (
                norm.cdf(t2c_full_ext[i + 1], s1mu, s1sd)
                - norm.cdf(t2c_full_ext[i], s1mu, s1sd)
            ) / c_area_rs1
            pi = (
                norm.cdf(t2c_full_ext[i + 1], s2mu, s2sd)
                - norm.cdf(t2c_full_ext[i], s2mu, s2sd)
            ) / i_area_rs1

            pc = max(pc, 1e-10)
            pi = max(pi, 1e-10)
            log_l += nc_rs1[i] * np.log(pc) + ni_rs1[i] * np.log(pi)

            pc2 = (
                (1 - norm.cdf(t2c_full_ext[n_ratings + i], s2mu, s2sd))
                - (1 - norm.cdf(t2c_full_ext[n_ratings + i + 1], s2mu, s2sd))
            ) / c_area_rs2
            pi2 = (
                (1 - norm.cdf(t2c_full_ext[n_ratings + i], s1mu, s1sd))
                - (1 - norm.cdf(t2c_full_ext[n_ratings + i + 1], s1mu, s1sd))
            ) / i_area_rs2

            pc2 = max(pc2, 1e-10)
            pi2 = max(pi2, 1e-10)
            log_l += nc_rs2[i] * np.log(pc2) + ni_rs2[i] * np.log(pi2)

        return -log_l

    opt = (
        {"maxiter": 10**5, "ftol": 1e-12}
        if matlab_compat
        else {"maxiter": 1000}
    )
    res = minimize(
        neg_log_likelihood,
        guess,
        method="SLSQP",
        bounds=bounds,
        constraints=cons,
        options=opt,
    )
    meta_d1_final = res.x[0]
    scale = np.sqrt(2 / (1 + s**2)) * s
    da = scale * d1
    meta_da = scale * meta_d1_final

    log_l = -res.fun
    # Free parameters: 1 (meta_d1) + (2*n_ratings - 2) type-2 criteria
    k = 2 * n_ratings - 1
    n_obs = int(np.sum(nr_s1) + np.sum(nr_s2))
    aic  = 2 * k - 2 * log_l
    bic  = k * np.log(n_obs) - 2 * log_l
    # AICc: small-sample correction
    aicc = aic + (2 * k * (k + 1)) / max(n_obs - k - 1, 1)

    # MATLAB packages the returned point even when fmincon quality is poor.
    success = True if matlab_compat else bool(res.success)

    return {
        "da": da,
        "meta_da": meta_da,
        "M_ratio": meta_da / da if da != 0 else np.nan,
        "M_diff": meta_da - da,
        "s": s,
        "logL": log_l,
        "k": k,
        "n": n_obs,
        "AIC": aic,
        "BIC": bic,
        "AICc": aicc,
        "success": success,
        "optimizer_success": bool(res.success),
    }

Type-2 measures

compute_type2_auc(nr_s1, nr_s2)

Compute Type-2 AUC from counts.

Source code in src/metasignal/stdpy/type2.py 
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
def compute_type2_auc(nr_s1: np.ndarray, nr_s2: np.ndarray) -> float:
    """Compute Type-2 AUC from counts."""
    n_ratings = len(nr_s1) // 2

    # correct/incorrect counts
    counts_c = nr_s2[n_ratings:] + nr_s1[:n_ratings][::-1]
    counts_i = nr_s1[n_ratings:] + nr_s2[:n_ratings][::-1]

    # Type-2 HR and FAR
    sum_c = np.sum(counts_c)
    sum_i = np.sum(counts_i)
    if sum_c == 0 or sum_i == 0:
        return np.nan
    hr2 = np.cumsum(counts_c[::-1]) / sum_c
    far2 = np.cumsum(counts_i[::-1]) / sum_i

    # Add (0,0) and sort
    hr2 = np.concatenate([[0], hr2])
    far2 = np.concatenate([[0], far2])

    # Trapezoidal integration
    auc = 0.0
    for i in range(len(hr2) - 1):
        auc += (far2[i + 1] - far2[i]) * (hr2[i + 1] + hr2[i]) / 2
    return auc

compute_gamma(nr_s1, nr_s2)

Compute Goodman-Kruskal's gamma.

Source code in src/metasignal/stdpy/type2.py 
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
def compute_gamma(nr_s1: np.ndarray, nr_s2: np.ndarray) -> float:
    """Compute Goodman-Kruskal's gamma."""
    n_ratings = len(nr_s1) // 2
    counts_i = nr_s1[n_ratings:] + nr_s2[:n_ratings][::-1]
    counts_c = nr_s2[n_ratings:] + nr_s1[:n_ratings][::-1]

    table = np.vstack([counts_i, counts_c])

    concordant = 0.0
    discordant = 0.0
    for i in range(2):
        for j in range(n_ratings):
            if i == 0:
                concordant += table[i, j] * np.sum(table[1, j + 1 :])
                discordant += table[i, j] * np.sum(table[1, :j])
            else:
                concordant += table[i, j] * np.sum(table[0, :j])
                discordant += table[i, j] * np.sum(table[0, j + 1 :])

    if (concordant + discordant) < 1e-10:
        return np.nan
    return (concordant - discordant) / (concordant + discordant)

compute_phi(nr_s1, nr_s2)

Compute Phi coefficient (correlation between correctness and confidence).

Source code in src/metasignal/stdpy/type2.py 
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
def compute_phi(nr_s1: np.ndarray, nr_s2: np.ndarray) -> float:
    """Compute Phi coefficient (correlation between correctness and confidence)."""
    n_ratings = len(nr_s1) // 2
    multiplier = np.concatenate([np.arange(n_ratings, 0, -1), np.arange(1, n_ratings + 1)])
    correct_trials_s1 = np.concatenate([np.ones(n_ratings), np.zeros(n_ratings)])
    correct_trials_s2 = np.concatenate([np.zeros(n_ratings), np.ones(n_ratings)])

    correct_cells = np.concatenate([nr_s1[:n_ratings], nr_s2[n_ratings:]])
    incorrect_cells = np.concatenate([nr_s2[:n_ratings], nr_s1[n_ratings:]])

    total = np.sum(correct_cells) + np.sum(incorrect_cells)
    if total == 0:
        return np.nan

    av_acc = np.sum(correct_cells) / total
    av_conf = np.sum((correct_cells + incorrect_cells) * multiplier) / total

    numerator = np.sum(
        (multiplier - av_conf) * (correct_trials_s1 - av_acc) * nr_s1 +
        (multiplier - av_conf) * (correct_trials_s2 - av_acc) * nr_s2
    )

    den1 = np.sum((multiplier - av_conf)**2 * (nr_s1 + nr_s2))
    den2 = np.sum((correct_trials_s1 - av_acc)**2 * (nr_s1 + nr_s2[::-1]))
    denominator = np.sqrt(den1 * den2)

    if denominator == 0:
        return np.nan

    return float(numerator / denominator)

compute_delta_conf(nr_s1, nr_s2)

Compute Delta Confidence and related expected measures.

Parameters:

Name Type Description Default
nr_s1 ndarray

Counts array for S1 stimulus responses

 required
nr_s2 ndarray

Counts array for S2 stimulus responses

 required

Returns:

Type Description
dict[str, float]

Dictionary containing delta_conf, delta_conf_ratio, and delta_conf_diff.
Source code in src/metasignal/stdpy/type2.py 
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
def compute_delta_conf(nr_s1: np.ndarray, nr_s2: np.ndarray) -> dict[str, float]:
    """Compute Delta Confidence and related expected measures.

    Args:
        nr_s1: Counts array for S1 stimulus responses
        nr_s2: Counts array for S2 stimulus responses

    Returns:
        Dictionary containing delta_conf, delta_conf_ratio, and delta_conf_diff.
    """
    n_ratings = len(nr_s1) // 2

    def _compute_delta(n_s1: np.ndarray, n_s2: np.ndarray) -> float:
        multiplier = np.concatenate(
            [np.arange(n_ratings, 0, -1), np.arange(1, n_ratings + 1)]
        )
        correct_cells = np.concatenate([n_s1[:n_ratings], n_s2[n_ratings:]])
        incorrect_cells = np.concatenate([n_s2[:n_ratings], n_s1[n_ratings:]])

        sum_c = np.sum(correct_cells)
        sum_i = np.sum(incorrect_cells)

        mean_conf_correct = np.sum(correct_cells * multiplier) / sum_c if sum_c > 0 else 0.0
        mean_conf_incorrect = np.sum(incorrect_cells * multiplier) / sum_i if sum_i > 0 else 0.0

        return float(mean_conf_correct - mean_conf_incorrect)

    # Actual delta conf
    delta_conf = _compute_delta(nr_s1, nr_s2)

    # Expected delta conf from SDT expectations
    sdt_expect = sdt_expect_conf(nr_s1, nr_s2)
    conf_diff_expected = _compute_delta(
        np.array(sdt_expect["nR_S1_exp"]),
        np.array(sdt_expect["nR_S2_exp"])
    )

    return {
        "delta_conf": delta_conf,
        "delta_conf_ratio": delta_conf / conf_diff_expected if conf_diff_expected != 0 else np.nan,
        "delta_conf_diff": delta_conf - conf_diff_expected,
    }

sdt_expect_conf(nr_s1, nr_s2)

Compute SDT-expected rating distributions (MATLAB SDTexpectConf).

Returns proportions for nR_S1_exp / nR_S2_exp (each sums to 1), matching MATLAB nR_*_SDTexpect. Using counts instead of proportions changes Ratio/Diff measures whenever S1 and S2 base rates differ (e.g. Locke).

Source code in src/metasignal/stdpy/type2.py 
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
def sdt_expect_conf(nr_s1: np.ndarray, nr_s2: np.ndarray) -> dict[str, Any]:
    """Compute SDT-expected rating distributions (MATLAB ``SDTexpectConf``).

    Returns **proportions** for ``nR_S1_exp`` / ``nR_S2_exp`` (each sums to 1),
    matching MATLAB ``nR_*_SDTexpect``. Using counts instead of proportions
    changes Ratio/Diff measures whenever S1 and S2 base rates differ (e.g. Locke).
    """
    n_ratings = len(nr_s1) // 2

    # Correct for empty cells
    if np.any(nr_s1 == 0) or np.any(nr_s2 == 0):
        nr_s1_c = nr_s1 + 1 / (2 * n_ratings)
        nr_s2_c = nr_s2 + 1 / (2 * n_ratings)
    else:
        nr_s1_c = nr_s1
        nr_s2_c = nr_s2

    sum1 = np.sum(nr_s1_c)
    sum2 = np.sum(nr_s2_c)

    # Match MATLAB: flip(cumsum(flip(nR(2:end)))) / sum(nR)
    # Excludes the first cell from the cumulative chain (same as SDTexpectConf.m).
    hr = np.flip(np.cumsum(np.flip(nr_s2_c[1:]))) / sum2
    far = np.flip(np.cumsum(np.flip(nr_s1_c[1:]))) / sum1

    # d' and c at the primary (middle) criterion — MATLAB HR(nRatings)
    t1_idx = n_ratings - 1
    dprime = norm.ppf(hr[t1_idx]) - norm.ppf(far[t1_idx])

    # c for all criteria
    c = -0.5 * (norm.ppf(hr) + norm.ppf(far))

    # Expected HR and FAR based on d' and c
    exp_far = 1 - norm.cdf(c, -dprime / 2, 1)
    exp_hr = 1 - norm.cdf(c, dprime / 2, 1)

    # Expected *proportions* (MATLAB nR_*_SDTexpect); do NOT rescale by trial counts
    exp_nr_s1 = np.flip(np.diff(np.concatenate([[0], np.flip(exp_far), [1]])))
    exp_nr_s2 = np.flip(np.diff(np.concatenate([[0], np.flip(exp_hr), [1]])))

    return {
        "nR_S1_exp": exp_nr_s1,
        "nR_S2_exp": exp_nr_s2,
        "nR_S1_act": nr_s1,
        "nR_S2_act": nr_s2,
        "dprime": dprime,
    }

Noise and uncertainty

compute_meta_noise(stim, resp, conf, n_ratings)

Fit Lognormal Meta Noise model (MATLAB-parity search).

Source code in src/metasignal/stdpy/metanoise.py 
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
def compute_meta_noise(
    stim: np.ndarray, resp: np.ndarray, conf: np.ndarray, n_ratings: int
) -> dict[str, Any]:
    """Fit Lognormal Meta Noise model (MATLAB-parity search)."""
    stim = np.asarray(stim)
    resp = np.asarray(resp)
    conf = np.asarray(conf)

    min_stim = int(np.min(stim))
    max_stim = int(np.max(stim))
    if min_stim == max_stim:
        raise ValueError("stim must contain both classes; found only one class in input.")

    nr_stim_s1 = np.zeros(2 * n_ratings)
    nr_stim_s2 = np.zeros(2 * n_ratings)
    for rating in range(1, n_ratings + 1):
        nr_stim_s1[n_ratings - rating] = np.sum(
            (stim == min_stim) & (resp == min_stim) & (conf == rating)
        )
        nr_stim_s2[n_ratings - rating] = np.sum(
            (stim == max_stim) & (resp == min_stim) & (conf == rating)
        )
        nr_stim_s1[n_ratings + rating - 1] = np.sum(
            (stim == min_stim) & (resp == max_stim) & (conf == rating)
        )
        nr_stim_s2[n_ratings + rating - 1] = np.sum(
            (stim == max_stim) & (resp == max_stim) & (conf == rating)
        )
    data_counts = np.vstack([nr_stim_s1, nr_stim_s2])

    dprime, c_init = _compute_sdt_criteria(stim, resp, conf, n_ratings)
    c_init = np.array(c_init, dtype=float)
    c_init[np.isnan(c_init)] = np.inf
    c = np.copy(c_init)

    # Refine each criterion on a local grid (matches MATLAB).
    for crit in range(2 * n_ratings - 1):
        if np.isinf(c_init[crit]):
            r_range = np.arange(-5.0, 5.0 + 1e-12, 0.01)
        else:
            r_range = np.arange(c_init[crit] - 0.5, c_init[crit] + 0.5 + 1e-12, 0.01)
        pr_s1 = norm.cdf(r_range, -dprime / 2.0, 1.0)
        pr_s2 = norm.cdf(r_range, dprime / 2.0, 1.0)
        pr_model = np.vstack([pr_s1, pr_s2, 1.0 - pr_s1, 1.0 - pr_s2])
        pr_model[pr_model == 0.0] = 1e-5
        data = np.concatenate(
            [
                np.sum(data_counts[:, : crit + 1], axis=1),
                np.sum(data_counts[:, crit + 1 :], axis=1),
            ]
        )
        logl_grid = -np.sum(np.log(pr_model) * data[:, None], axis=0)
        c[crit] = r_range[int(np.argmin(logl_grid))]

    # Gaussian baseline (metaNoise = 0).
    # MATLAB uses diff([0, normcdf(c), 1]) — equivalent for finite c.
    pr_s1 = np.diff(np.concatenate([[0.0], norm.cdf(c, -dprime / 2.0, 1.0), [1.0]]))
    pr_s2 = np.diff(np.concatenate([[0.0], norm.cdf(c, dprime / 2.0, 1.0), [1.0]]))
    pr_model = np.vstack([pr_s1, pr_s2])
    prob_dec_resp = np.column_stack(
        [
            np.sum(pr_model[:, :n_ratings], axis=1),
            np.sum(pr_model[:, n_ratings:], axis=1),
        ]
    )
    pr_clip = pr_model.copy()
    pr_clip[pr_clip == 0.0] = 1e-5
    logl_mn0 = float(-np.sum(np.log(pr_clip) * data_counts))

    # Seed criterion cache at metaNoise=0 (log-space of Type-1 criteria).
    pos_crit = c[n_ratings:] - c[n_ratings - 1]
    neg_crit = -(c[n_ratings - 2 :: -1] - c[n_ratings - 1])
    crit_values_all = np.vstack([pos_crit, neg_crit])
    # Tiny FP negatives from criterion-grid differencing are treated as zero
    # (same convention as the pr_full clipping below; avoids log of a negative).
    crit_values_all[crit_values_all < 0.0] = 0.0
    crit_values_all[crit_values_all == 0.0] = 1e-5
    crit_lognorm = np.log(crit_values_all)

    state = _FitState(
        meta_noise_tested=[0.0],
        criteria_for_tested=[crit_lognorm],
    )

    def _mn_func(mn: float) -> dict[str, Any]:
        return _logl_func_meta_noise(
            float(mn), data_counts, prob_dec_resp, dprime, c, state
        )

    lower_bound_info = {"x": 0.0, "logL": logl_mn0}
    best = _search_with_lower_bound(lower_bound_info, _mn_func, is_meta_noise=True)

    return {
        "meta_noise": float(best["x"]),
        "dprime": float(dprime),
        "c": c.tolist(),
        "logL": float(best["logL"]),
    }

compute_meta_uncertainty(stim, resp, conf, n_ratings, rng=None, n_starts=5, *, matlab_compat=False)

Estimate meta-uncertainty (MATLAB compute_metaUncertainty).

Parameters

rng : Optional generator. Used for an extra start in multi-start mode, or as the sole MATLAB-style start when matlab_compat=True. n_starts : Number of deterministic optimizer starts (default 5). Ignored when matlab_compat=True. matlab_compat : If True, mimic MATLAB's single random start [1, 0, 0.2, sort(2*rand(...))] (requires rng or uses a fresh default generator). Prefer multi-start for stable science; use this only when reproducing a specific MATLAB Result file.

Source code in src/metasignal/stdpy/uncertainty.py 
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
def compute_meta_uncertainty(
    stim: np.ndarray,
    resp: np.ndarray,
    conf: np.ndarray,
    n_ratings: int,
    rng: np.random.Generator | None = None,
    n_starts: int = 5,
    *,
    matlab_compat: bool = False,
) -> float:
    """Estimate meta-uncertainty (MATLAB compute_metaUncertainty).

    Parameters
    ----------
    rng :
        Optional generator. Used for an extra start in multi-start mode, or
        as the sole MATLAB-style start when ``matlab_compat=True``.
    n_starts :
        Number of deterministic optimizer starts (default 5). Ignored when
        ``matlab_compat=True``.
    matlab_compat :
        If True, mimic MATLAB's single random start
        ``[1, 0, 0.2, sort(2*rand(...))]`` (requires ``rng`` or uses a
        fresh default generator). Prefer multi-start for stable science;
        use this only when reproducing a specific MATLAB Result file.
    """
    stim_values = np.unique(stim)
    n_stim = len(stim_values)
    if n_stim < 2:
        raise ValueError("stim must contain both classes; found only one class in input.")

    n_choice = np.zeros((n_ratings * 2, n_stim))
    resp_min = np.min(resp)
    resp_max = np.max(resp)

    for i, sv in enumerate(stim_values):
        mask = stim == sv
        for r in range(1, n_ratings + 1):
            n_choice[n_ratings - r, i] = np.sum(
                (resp[mask] == resp_min) & (conf[mask] == r)
            )
        for r in range(1, n_ratings + 1):
            n_choice[n_ratings + r - 1, i] = np.sum(
                (resp[mask] == resp_max) & (conf[mask] == r)
            )

    def objective(params: np.ndarray) -> float:
        llh = get_llh_choice(stim_values, params)
        # MATLAB does not clip; guess_rate keeps llh > 0 in practice.
        if matlab_compat:
            llh = np.maximum(llh, np.finfo(float).tiny)
        else:
            llh = np.clip(llh, 1e-10, 1.0)
        return float(-np.sum(n_choice * np.log(llh)))

    n_conf_crit = n_ratings - 1
    bounds = [
        (0.0, 10.0),  # stim_sens
        (-3.0, 3.0),  # stim_crit
        (0.01, 5.0),  # meta_uncertainty
    ] + [(0.0, 5.0)] * n_conf_crit

    if matlab_compat:
        rng_m = rng if rng is not None else np.random.default_rng()
        starts = [
            np.concatenate([[1.0, 0.0, 0.2], np.sort(2.0 * rng_m.random(n_conf_crit))])
        ]
        opt = {"maxiter": 10**5, "maxfun": 10**5, "ftol": 1e-12}
    else:
        starts = _starting_points(n_conf_crit, max(1, int(n_starts)))
        if rng is not None:
            starts.append(
                np.concatenate(
                    [[1.0, 0.0, 0.2], np.sort(2.0 * rng.random(n_conf_crit))]
                )
            )
        opt = {"maxiter": 2000, "maxfun": 4000, "ftol": 1e-10}

    best_val = np.inf
    best_x = starts[0].copy()
    for guess in starts:
        res = minimize(
            objective,
            guess,
            method="L-BFGS-B",
            bounds=bounds,
            options=opt,
        )
        if np.isfinite(res.fun) and res.fun < best_val:
            best_val = float(res.fun)
            best_x = np.asarray(res.x, dtype=float)

    return float(best_x[2])

Group-level fitting

fit_group(data=None, stimuli='Stimuli', responses='Responses', confidence='Confidence', nRatings=4, subject=None, within=None, between=None, measures=None, nR_S1=None, nR_S2=None, method='mle', **bayes_kwargs)

Compute type-1 SDT and type-2 metacognitive measures for every group cell.

Parameters

data : pd.DataFrame Trial-level DataFrame with columns for stimuli, responses, and confidence. Subject / condition columns are optional. stimuli : str Column for stimulus identity (0 = S1, 1 = S2). responses : str Column for participant response (0 = S1, 1 = S2). confidence : str Column for confidence rating (integer, 1..nRatings). nRatings : int Number of confidence rating levels. subject : str or None Column identifying participants. None treats all rows as one cell. within : str or list of str or None Within-subjects factor column(s). When a single within-subjects factor with exactly two levels is supplied together with method='bayesian', fit_within_subject_comparison is used automatically. between : str or list of str or None Between-subjects factor column(s). When supplied together with method='bayesian', fit_full_metad_comparison is used. nR_S1 : array-like or None Count vector for S1 stimulus, length 2 * nRatings. When provided with nR_S2, data is ignored and a single-cell fit is run directly from counts. Equivalent to metadpy's metad(nR_S1=..., nR_S2=...). nR_S2 : array-like or None Count vector for S2 stimulus (same layout as nR_S1). method : str Fitting method. One of:

``'mle'`` *(default)*
    Frequentist MLE via :func:`fit_meta_d_mle`.  Returns a tidy
    DataFrame with one row per cell and 28 measure columns.

``'bayesian'``
    Hierarchical Bayesian via ``metasignal.sdtbayes``.  Requires
    ``pip install "metasignal[sdtbayes] @ git+https://github.com/saurabhr/metasignal.git"``
    and ``setup_runtime()``.
    The specific sdtbayes function is chosen automatically:

    - No ``subject`` → ``fit_subject_level`` (single participant)
    - ``subject`` only → ``fit_full_metad`` (group, no comparison)
    - ``between`` set → ``fit_full_metad_comparison`` (between-groups)
    - ``within`` with 2 levels → ``fit_within_subject_comparison``
    - ``within`` with >2 levels → ``fit_hierarchical_metad``

    Returns the raw ``FitResult`` from sdtbayes (use
    ``metasignal.sdtbayes.posterior_summary`` to inspect it).

``'two_stage'``
    Fast two-stage Bayesian: MLE per participant then Bayesian
    pooling via ``fit_two_stage_group`` / ``fit_two_stage_comparison``.

``'robust'``
    Full HMeta-d with Student-t hyperprior (downweights outliers).

``'hierarchical'``
    Ordered-logistic brms model (trial-level, extensible with
    ``items=`` for crossed effects).

**bayes_kwargs Extra keyword arguments forwarded to the chosen sdtbayes function (e.g. draws=2000, chains=4, target_accept=0.9, items=[...]).

Returns

results : pd.DataFrame or FitResult When method='mle': tidy DataFrame, one row per group cell, 28 measure columns. When any Bayesian method is chosen: the raw FitResult from metasignal.sdtbayes (use metasignal.sdtbayes.posterior_summary to inspect posteriors).

Examples

MLE (default), all measures, single subject:

from metasignal.stdpy import trialSimulation, fit_group df = trialSimulation(d=1.5, metad=1.5, nTrials=300) fit_group(df, nRatings=4).round(3)

Subject-level Bayesian (metadpy hmetad equivalent):

fit_group(nR_S1=nR_S1, nR_S2=nR_S2, nRatings=4, method='bayesian')

Group Bayesian — between-groups:

fit_group(df, subject='Subject', between='Group', ... nRatings=4, method='bayesian', draws=2000)

Within-subject Bayesian comparison:

fit_group(df, subject='Subject', within='Condition', ... nRatings=4, method='bayesian')

Only type-1 measures across subjects:

from metasignal.stdpy import responseSimulation, fit_group df = responseSimulation(d=1.5, metad=1.5, nSubjects=10, nTrials=200) fit_group(df, subject='Subject', nRatings=4, ... measures=['dprime', 'criterion', 'hit_rate', 'fa_rate'])

metadpy-style subset (meta-d', M_ratio, M_diff) with conditions:

from metasignal.stdpy import pairedResponseSimulation, fit_group df = pairedResponseSimulation(nSubjects=10, nTrials=200) fit_group(df, subject='Subject', within='Condition', nRatings=4, ... measures=['dprime', 'meta_d', 'M_ratio', 'M_diff'])

Source code in src/metasignal/stdpy/group.py 
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
def fit_group(
    data: pd.DataFrame = None,
    stimuli: str = "Stimuli",
    responses: str = "Responses",
    confidence: str = "Confidence",
    nRatings: int = 4,
    subject: Optional[str] = None,
    within: Optional[Union[str, List[str]]] = None,
    between: Optional[Union[str, List[str]]] = None,
    measures: Optional[Union[str, List[str]]] = None,
    nR_S1: Optional[np.ndarray] = None,
    nR_S2: Optional[np.ndarray] = None,
    method: str = "mle",
    **bayes_kwargs,
) -> "pd.DataFrame | object":
    """Compute type-1 SDT and type-2 metacognitive measures for every group cell.

    Parameters
    ----------
    data : pd.DataFrame
        Trial-level DataFrame with columns for stimuli, responses, and
        confidence. Subject / condition columns are optional.
    stimuli : str
        Column for stimulus identity (0 = S1, 1 = S2).
    responses : str
        Column for participant response (0 = S1, 1 = S2).
    confidence : str
        Column for confidence rating (integer, 1..nRatings).
    nRatings : int
        Number of confidence rating levels.
    subject : str or None
        Column identifying participants. ``None`` treats all rows as one cell.
    within : str or list of str or None
        Within-subjects factor column(s).  When a single within-subjects
        factor with exactly two levels is supplied together with
        ``method='bayesian'``, ``fit_within_subject_comparison`` is used
        automatically.
    between : str or list of str or None
        Between-subjects factor column(s).  When supplied together with
        ``method='bayesian'``, ``fit_full_metad_comparison`` is used.
    nR_S1 : array-like or None
        Count vector for S1 stimulus, length ``2 * nRatings``.
        When provided with ``nR_S2``, ``data`` is ignored and a single-cell
        fit is run directly from counts.
        Equivalent to metadpy's ``metad(nR_S1=..., nR_S2=...)``.
    nR_S2 : array-like or None
        Count vector for S2 stimulus (same layout as ``nR_S1``).
    method : str
        Fitting method.  One of:

        ``'mle'`` *(default)*
            Frequentist MLE via :func:`fit_meta_d_mle`.  Returns a tidy
            DataFrame with one row per cell and 28 measure columns.

        ``'bayesian'``
            Hierarchical Bayesian via ``metasignal.sdtbayes``.  Requires
            ``pip install "metasignal[sdtbayes] @ git+https://github.com/saurabhr/metasignal.git"``
            and ``setup_runtime()``.
            The specific sdtbayes function is chosen automatically:

            - No ``subject`` → ``fit_subject_level`` (single participant)
            - ``subject`` only → ``fit_full_metad`` (group, no comparison)
            - ``between`` set → ``fit_full_metad_comparison`` (between-groups)
            - ``within`` with 2 levels → ``fit_within_subject_comparison``
            - ``within`` with >2 levels → ``fit_hierarchical_metad``

            Returns the raw ``FitResult`` from sdtbayes (use
            ``metasignal.sdtbayes.posterior_summary`` to inspect it).

        ``'two_stage'``
            Fast two-stage Bayesian: MLE per participant then Bayesian
            pooling via ``fit_two_stage_group`` / ``fit_two_stage_comparison``.

        ``'robust'``
            Full HMeta-d with Student-t hyperprior (downweights outliers).

        ``'hierarchical'``
            Ordered-logistic brms model (trial-level, extensible with
            ``items=`` for crossed effects).

    **bayes_kwargs
        Extra keyword arguments forwarded to the chosen sdtbayes function
        (e.g. ``draws=2000``, ``chains=4``, ``target_accept=0.9``,
        ``items=[...]``).

    Returns
    -------
    results : pd.DataFrame or FitResult
        When ``method='mle'``: tidy DataFrame, one row per group cell,
        28 measure columns.  When any Bayesian method is chosen: the raw
        ``FitResult`` from ``metasignal.sdtbayes`` (use
        ``metasignal.sdtbayes.posterior_summary`` to inspect posteriors).

    Examples
    --------
    MLE (default), all measures, single subject:

    >>> from metasignal.stdpy import trialSimulation, fit_group
    >>> df = trialSimulation(d=1.5, metad=1.5, nTrials=300)
    >>> fit_group(df, nRatings=4).round(3)

    Subject-level Bayesian (metadpy ``hmetad`` equivalent):

    >>> fit_group(nR_S1=nR_S1, nR_S2=nR_S2, nRatings=4, method='bayesian')

    Group Bayesian — between-groups:

    >>> fit_group(df, subject='Subject', between='Group',
    ...           nRatings=4, method='bayesian', draws=2000)

    Within-subject Bayesian comparison:

    >>> fit_group(df, subject='Subject', within='Condition',
    ...           nRatings=4, method='bayesian')

    Only type-1 measures across subjects:

    >>> from metasignal.stdpy import responseSimulation, fit_group
    >>> df = responseSimulation(d=1.5, metad=1.5, nSubjects=10, nTrials=200)
    >>> fit_group(df, subject='Subject', nRatings=4,
    ...           measures=['dprime', 'criterion', 'hit_rate', 'fa_rate'])

    metadpy-style subset (meta-d', M_ratio, M_diff) with conditions:

    >>> from metasignal.stdpy import pairedResponseSimulation, fit_group
    >>> df = pairedResponseSimulation(nSubjects=10, nTrials=200)
    >>> fit_group(df, subject='Subject', within='Condition', nRatings=4,
    ...           measures=['dprime', 'meta_d', 'M_ratio', 'M_diff'])
    """
    keep_cols = _resolve_measures(measures)

    # --- Bayesian routing ---
    if method != "mle":
        return _fit_group_bayesian(
            data=data, stimuli=stimuli, responses=responses,
            confidence=confidence, nRatings=nRatings,
            subject=subject, within=within, between=between,
            nR_S1=nR_S1, nR_S2=nR_S2,
            method=method, **bayes_kwargs,
        )

    # --- Count-array shortcut (single cell, MLE, no grouping) ---
    if nR_S1 is not None and nR_S2 is not None:
        from metasignal.stdpy.simulate import ratings2df
        nr1 = np.asarray(nR_S1, dtype=float)
        nr2 = np.asarray(nR_S2, dtype=float)
        if len(nr1) != 2 * nRatings:
            raise ValueError(f"nR_S1 length {len(nr1)} != 2 * nRatings ({2 * nRatings})")
        df_counts = ratings2df(nr1, nr2)
        row = _fit_cell(df_counts, "Stimuli", "Responses", "Confidence", nRatings)
        return pd.DataFrame([{k: v for k, v in row.items() if k in keep_cols}])

    if data is None:
        raise ValueError("Provide either 'data' (DataFrame) or both 'nR_S1' and 'nR_S2'.")

    within_cols:  List[str] = _to_list(within)
    between_cols: List[str] = _to_list(between)

    group_cols: List[str] = []
    if subject is not None:
        group_cols.append(subject)
    group_cols.extend(within_cols)
    group_cols.extend(between_cols)

    if not group_cols:
        row = _fit_cell(data, stimuli, responses, confidence, nRatings)
        return pd.DataFrame([{k: v for k, v in row.items() if k in keep_cols}])

    rows = []
    for keys, grp in data.groupby(group_cols, sort=True):
        if not isinstance(keys, tuple):
            keys = (keys,)
        group_meta = dict(zip(group_cols, keys))
        cell = _fit_cell(grp, stimuli, responses, confidence, nRatings)
        rows.append({**group_meta, **{k: v for k, v in cell.items() if k in keep_cols}})

    return pd.DataFrame(rows).reset_index(drop=True)

Plotting

Note: stdpy.plot_forest and sdtbayes.plot_forest are unrelated functions that happen to share a name — the former plots any scalar measure across subjects/conditions, the latter plots MCMC posterior forest plots.

plot_confidence(nR_S1, nR_S2, ax=None, title='Confidence rating distribution', xlabel='Rating (S1 ← → S2)')

Bar chart of response-count distributions split by stimulus class.

Parameters

nR_S1, nR_S2 : Count vectors of length 2 * nRatings. ax : Existing axes to draw on. If None a new figure is created. title, xlabel : Axis labels.

Returns

ax : plt.Axes

Source code in src/metasignal/stdpy/plot.py 
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
def plot_confidence(
    nR_S1: np.ndarray,
    nR_S2: np.ndarray,
    ax: Optional[plt.Axes] = None,
    title: str = "Confidence rating distribution",
    xlabel: str = "Rating (S1 ← → S2)",
) -> plt.Axes:
    """Bar chart of response-count distributions split by stimulus class.

    Parameters
    ----------
    nR_S1, nR_S2 :
        Count vectors of length ``2 * nRatings``.
    ax :
        Existing axes to draw on. If ``None`` a new figure is created.
    title, xlabel :
        Axis labels.

    Returns
    -------
    ax : plt.Axes
    """
    nR_S1 = np.asarray(nR_S1, dtype=float)
    nR_S2 = np.asarray(nR_S2, dtype=float)
    nRatings = len(nR_S1) // 2

    if ax is None:
        _, ax = plt.subplots(figsize=(8, 4))

    x = np.arange(2 * nRatings)
    rating_labels = (
        [f"S1-r{i}" for i in range(nRatings, 0, -1)]
        + [f"S2-r{i}" for i in range(1, nRatings + 1)]
    )

    # S1 stim: correct (S1 response half) and incorrect (S2 response half)
    ax.bar(x[:nRatings],  nR_S1[:nRatings],  color=_S1_COLOR, alpha=_CORR_ALPHA,  label="S1 stim, S1 resp (correct)")
    ax.bar(x[nRatings:],  nR_S1[nRatings:],  color=_S1_COLOR, alpha=_INCORR_ALPHA, label="S1 stim, S2 resp (error)")
    # S2 stim: incorrect (S1 response half) and correct (S2 response half)
    ax.bar(x[:nRatings],  nR_S2[:nRatings],  color=_S2_COLOR, alpha=_INCORR_ALPHA, label="S2 stim, S1 resp (error)",
           bottom=nR_S1[:nRatings])
    ax.bar(x[nRatings:],  nR_S2[nRatings:],  color=_S2_COLOR, alpha=_CORR_ALPHA,  label="S2 stim, S2 resp (correct)",
           bottom=nR_S1[nRatings:])

    ax.axvline(nRatings - 0.5, color="black", linewidth=1.2, linestyle="--")
    ax.set_xticks(x)
    ax.set_xticklabels(rating_labels, rotation=45, ha="right", fontsize=8)
    ax.set_ylabel("Count")
    ax.set_xlabel(xlabel)
    ax.set_title(title)
    ax.legend(fontsize=8, loc="upper center", ncol=2)
    return ax

plot_type2roc(nR_S1, nR_S2, ax=None, show_ideal=True, label='Observed', color=_S2_COLOR, title='Type-2 ROC')

Plot the Type-2 ROC curve.

Parameters

nR_S1, nR_S2 : Count vectors of length 2 * nRatings. ax : Existing axes. None creates a new figure. show_ideal : Overlay the ideal SDT-predicted ROC curve. label : Legend label for the observed curve. color : Line / marker colour. title : Axes title.

Returns

ax : plt.Axes

Source code in src/metasignal/stdpy/plot.py 
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
def plot_type2roc(
    nR_S1: np.ndarray,
    nR_S2: np.ndarray,
    ax: Optional[plt.Axes] = None,
    show_ideal: bool = True,
    label: str = "Observed",
    color: str = _S2_COLOR,
    title: str = "Type-2 ROC",
) -> plt.Axes:
    """Plot the Type-2 ROC curve.

    Parameters
    ----------
    nR_S1, nR_S2 :
        Count vectors of length ``2 * nRatings``.
    ax :
        Existing axes. ``None`` creates a new figure.
    show_ideal :
        Overlay the ideal SDT-predicted ROC curve.
    label :
        Legend label for the observed curve.
    color :
        Line / marker colour.
    title :
        Axes title.

    Returns
    -------
    ax : plt.Axes
    """
    nR_S1 = np.asarray(nR_S1, dtype=float)
    nR_S2 = np.asarray(nR_S2, dtype=float)
    nRatings = len(nR_S1) // 2

    if ax is None:
        _, ax = plt.subplots(figsize=(5, 5))

    # Observed type-2 ROC points
    counts_c = nR_S2[nRatings:] + nR_S1[:nRatings][::-1]
    counts_i = nR_S1[nRatings:] + nR_S2[:nRatings][::-1]
    hr2  = np.concatenate([[0], np.cumsum(counts_c[::-1]) / (counts_c.sum() or 1)])
    far2 = np.concatenate([[0], np.cumsum(counts_i[::-1]) / (counts_i.sum() or 1)])
    auc2 = compute_type2_auc(nR_S1, nR_S2)

    ax.plot(far2, hr2, "o-", color=color, label=f"{label} (AUC2={auc2:.3f})")

    # Ideal SDT overlay
    if show_ideal:
        sdt = sdt_expect_conf(nR_S1, nR_S2)
        nr1e = np.array(sdt["nR_S1_exp"])
        nr2e = np.array(sdt["nR_S2_exp"])
        c_e  = nr2e[nRatings:] + nr1e[:nRatings][::-1]
        i_e  = nr1e[nRatings:] + nr2e[:nRatings][::-1]
        hr2e  = np.concatenate([[0], np.cumsum(c_e[::-1]) / (c_e.sum() or 1)])
        far2e = np.concatenate([[0], np.cumsum(i_e[::-1]) / (i_e.sum() or 1)])
        auc2e = compute_type2_auc(nr1e, nr2e)
        ax.plot(far2e, hr2e, "s--", color="grey", alpha=0.7,
                label=f"SDT ideal (AUC2={auc2e:.3f})")

    ax.plot([0, 1], [0, 1], "k:", linewidth=1, label="Chance")
    ax.set_xlabel("Type-2 False Alarm Rate")
    ax.set_ylabel("Type-2 Hit Rate")
    ax.set_xlim(0, 1); ax.set_ylim(0, 1)
    ax.set_title(title)
    ax.legend(fontsize=8)
    ax.set_aspect("equal")
    return ax

plot_sanity_check(nR_S1, nR_S2, meta_d=None, title='Sanity check')

Three-panel diagnostic for a single participant.

Panels: (1) confidence distribution, (2) type-2 ROC, (3) text summary of key SDT parameters.

Parameters

nR_S1, nR_S2 : Count vectors of length 2 * nRatings. meta_d : Optional pre-computed meta-d' to include in the summary panel. title : Figure suptitle.

Returns

fig : plt.Figure axes : np.ndarray of plt.Axes (shape (3,))

Source code in src/metasignal/stdpy/plot.py 
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
def plot_sanity_check(
    nR_S1: np.ndarray,
    nR_S2: np.ndarray,
    meta_d: Optional[float] = None,
    title: str = "Sanity check",
) -> Tuple[plt.Figure, np.ndarray]:
    """Three-panel diagnostic for a single participant.

    Panels: (1) confidence distribution, (2) type-2 ROC,
    (3) text summary of key SDT parameters.

    Parameters
    ----------
    nR_S1, nR_S2 :
        Count vectors of length ``2 * nRatings``.
    meta_d :
        Optional pre-computed meta-d' to include in the summary panel.
    title :
        Figure suptitle.

    Returns
    -------
    fig : plt.Figure
    axes : np.ndarray of plt.Axes  (shape (3,))
    """
    nR_S1 = np.asarray(nR_S1, dtype=float)
    nR_S2 = np.asarray(nR_S2, dtype=float)
    nRatings = len(nR_S1) // 2

    fig, axes = plt.subplots(1, 3, figsize=(15, 4))

    # Panel 1 — confidence distribution
    plot_confidence(nR_S1, nR_S2, ax=axes[0], title="Confidence distribution")

    # Panel 2 — type-2 ROC
    plot_type2roc(nR_S1, nR_S2, ax=axes[1], title="Type-2 ROC")

    # Panel 3 — parameter summary text
    ax3 = axes[2]
    ax3.axis("off")
    sdt = sdt_expect_conf(nR_S1, nR_S2)
    dprime = sdt["dprime"]
    auc2   = compute_type2_auc(nR_S1, nR_S2)
    n_total = int(nR_S1.sum() + nR_S2.sum())

    lines = [
        f"N trials : {n_total}",
        f"d'       : {dprime:.3f}",
        f"AUC2     : {auc2:.3f}",
    ]
    if meta_d is not None:
        lines.append(f"meta-d'  : {meta_d:.3f}")
        lines.append(f"M_ratio  : {meta_d / dprime:.3f}" if dprime > 0 else "M_ratio  : N/A")

    text = "\n".join(lines)
    ax3.text(0.1, 0.55, text, transform=ax3.transAxes,
             fontsize=11, verticalalignment="center", fontfamily="monospace",
             bbox=dict(boxstyle="round", facecolor="whitesmoke", alpha=0.8))
    ax3.set_title("Summary")

    fig.suptitle(title, fontsize=13, y=1.02)
    fig.tight_layout()
    return fig, axes

plot_forest(data, measure, group_col=None, subject_col='Subject', ci=0.95, ax=None, title=None, palette=None)

Forest plot of a scalar measure across subjects / conditions.

Parameters

data : pd.DataFrame Output of :func:fit_group — one row per subject (× condition). measure : str Column name of the measure to plot (e.g. 'M_ratio'). group_col : str or None Column to colour-code groups / conditions. None = single colour. subject_col : str or None Column identifying individual subjects (plotted as points). ci : float Confidence interval coverage for the group mean bar (default 95 %). ax : Existing axes. None creates a new figure. title : str or None Axes title (defaults to measure name). palette : list of str or None Colours per group level.

Returns

ax : plt.Axes

Source code in src/metasignal/stdpy/plot.py 
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
def plot_forest(
    data,
    measure: str,
    group_col: Optional[str] = None,
    subject_col: Optional[str] = "Subject",
    ci: float = 0.95,
    ax: Optional[plt.Axes] = None,
    title: Optional[str] = None,
    palette: Optional[List[str]] = None,
) -> plt.Axes:
    """Forest plot of a scalar measure across subjects / conditions.

    Parameters
    ----------
    data : pd.DataFrame
        Output of :func:`fit_group` — one row per subject (× condition).
    measure : str
        Column name of the measure to plot (e.g. ``'M_ratio'``).
    group_col : str or None
        Column to colour-code groups / conditions. ``None`` = single colour.
    subject_col : str or None
        Column identifying individual subjects (plotted as points).
    ci : float
        Confidence interval coverage for the group mean bar (default 95 %).
    ax :
        Existing axes. ``None`` creates a new figure.
    title : str or None
        Axes title (defaults to measure name).
    palette : list of str or None
        Colours per group level.

    Returns
    -------
    ax : plt.Axes
    """
    import pandas as _pd

    data = _pd.DataFrame(data)
    groups = [None] if group_col is None else sorted(data[group_col].unique())
    palette = palette or [_S1_COLOR, _S2_COLOR, "#4CAF50", "#9C27B0"]

    if ax is None:
        _, ax = plt.subplots(figsize=(max(4, len(groups) * 2.5), 4))

    for gi, grp in enumerate(groups):
        color = palette[gi % len(palette)]
        sub = data if grp is None else data[data[group_col] == grp]
        vals = sub[measure].dropna().values
        if len(vals) == 0:
            continue

        mean = vals.mean()
        se   = vals.std(ddof=1) / np.sqrt(len(vals))
        tcrit = t_dist.ppf((1 + ci) / 2, df=len(vals) - 1) if len(vals) > 1 else 0
        ci_lo, ci_hi = mean - tcrit * se, mean + tcrit * se

        x = gi
        label = str(grp) if grp is not None else measure

        # Individual subject points
        jitter = np.random.default_rng(gi).uniform(-0.15, 0.15, len(vals))
        ax.scatter(np.full(len(vals), x) + jitter, vals,
                   color=color, alpha=0.4, s=25, zorder=2)

        # Mean ± CI
        ax.errorbar(x, mean, yerr=[[mean - ci_lo], [ci_hi - mean]],
                    fmt="D", color=color, markersize=8, linewidth=2,
                    capsize=5, zorder=3, label=label)

    ax.axhline(0, color="black", linewidth=0.8, linestyle=":")
    ax.set_xticks(range(len(groups)))
    ax.set_xticklabels([str(g) for g in groups] if groups[0] is not None else [""])
    ax.set_ylabel(measure)
    ax.set_title(title or measure)
    if group_col is not None:
        ax.legend(title=group_col, fontsize=8)
    return ax

plot_measures(data, measures, group_col=None, ci=0.95, palette=None, figsize=None)

Bar chart comparing one or more measures across conditions / groups.

Parameters

data : pd.DataFrame Output of :func:fit_group. measures : str or list of str Measure column(s) to plot. One subplot per measure. group_col : str or None Column used to split bars by colour. ci : float Confidence interval coverage (default 95 %). palette : list of str or None Colours per group level. figsize : tuple or None Figure size override.

Returns

fig : plt.Figure axes : np.ndarray of plt.Axes

Source code in src/metasignal/stdpy/plot.py 
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
def plot_measures(
    data,
    measures: Union[str, List[str]],
    group_col: Optional[str] = None,
    ci: float = 0.95,
    palette: Optional[List[str]] = None,
    figsize: Optional[Tuple[float, float]] = None,
) -> Tuple[plt.Figure, np.ndarray]:
    """Bar chart comparing one or more measures across conditions / groups.

    Parameters
    ----------
    data : pd.DataFrame
        Output of :func:`fit_group`.
    measures : str or list of str
        Measure column(s) to plot. One subplot per measure.
    group_col : str or None
        Column used to split bars by colour.
    ci : float
        Confidence interval coverage (default 95 %).
    palette : list of str or None
        Colours per group level.
    figsize : tuple or None
        Figure size override.

    Returns
    -------
    fig : plt.Figure
    axes : np.ndarray of plt.Axes
    """
    import pandas as _pd

    data = _pd.DataFrame(data)
    if isinstance(measures, str):
        measures = [measures]

    groups  = [None] if group_col is None else sorted(data[group_col].unique())
    palette = palette or [_S1_COLOR, _S2_COLOR, "#4CAF50", "#9C27B0"]
    n_meas  = len(measures)
    n_grp   = len(groups)

    fig, axes = plt.subplots(1, n_meas,
                             figsize=figsize or (max(4, n_meas * 3), 4),
                             sharey=False)
    if n_meas == 1:
        axes = np.array([axes])

    bar_width = 0.7 / max(n_grp, 1)

    for mi, meas in enumerate(measures):
        ax = axes[mi]
        for gi, grp in enumerate(groups):
            color = palette[gi % len(palette)]
            sub   = data if grp is None else data[data[group_col] == grp]
            vals  = sub[meas].dropna().values
            if len(vals) == 0:
                continue

            mean  = vals.mean()
            se    = vals.std(ddof=1) / np.sqrt(len(vals)) if len(vals) > 1 else 0
            tcrit = t_dist.ppf((1 + ci) / 2, df=max(len(vals) - 1, 1))
            err   = tcrit * se

            xpos = gi * bar_width - (n_grp - 1) * bar_width / 2
            label = str(grp) if grp is not None else meas
            ax.bar(xpos, mean, width=bar_width * 0.85,
                   color=color, alpha=0.8, label=label, zorder=2)
            ax.errorbar(xpos, mean, yerr=err,
                        fmt="none", color="black", linewidth=1.5, capsize=4, zorder=3)

        ax.axhline(0, color="black", linewidth=0.6, linestyle=":")
        ax.set_title(meas)
        ax.set_xticks([])
        if mi == 0:
            ax.set_ylabel("Value")
        if group_col is not None and mi == n_meas - 1:
            ax.legend(title=group_col, fontsize=8)

    fig.tight_layout()
    return fig, axes

Simulation

type2_SDT_simuation(d=1.0, noise=0.2, c=0.0, nRatings=4, nTrials=500, rng=None)

Simulate Type-2 SDT response counts.

Generates nR_S1 and nR_S2 count arrays compatible with :func:metasignal.stdpy.fit_meta_d_mle and :func:metasignal.stdpy.trials_to_counts.

Parameters

d : Type-1 d-prime (sensitivity). noise : Standard deviation of metacognitive noise added to the type-1 internal response for the type-2 judgment. Pass a two-element list [sigma_rS1, sigma_rS2] to simulate response-conditional noise. c : Type-1 decision criterion. nRatings : Number of confidence rating levels. nTrials : Number of trials to simulate. rng : Optional :class:numpy.random.Generator for reproducibility. If None, uses :func:numpy.random.default_rng.

Returns

nR_S1, nR_S2 : ndarray Response-count vectors of length 2 * nRatings.

*Ordering (same convention as metadpy / Maniscalco & Lau):*

``nR_S1[0 : nRatings]``   — S1 stimulus, S1 response, ratings high→low
``nR_S1[nRatings : 2*nRatings]`` — S1 stimulus, S2 response, ratings low→high

Same layout for nR_S2.

Examples

nR_S1, nR_S2 = type2_SDT_simuation(d=1.5, noise=0.0, c=0, nRatings=4, nTrials=1000) len(nR_S1) 8

Source code in src/metasignal/stdpy/simulate.py 
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
def type2_SDT_simuation(
    d: float = 1.0,
    noise: Union[float, list[float]] = 0.2,
    c: float = 0.0,
    nRatings: int = 4,
    nTrials: int = 500,
    rng: np.random.Generator | None = None,
) -> tuple[np.ndarray, np.ndarray]:
    """Simulate Type-2 SDT response counts.

    Generates nR_S1 and nR_S2 count arrays compatible with
    :func:`metasignal.stdpy.fit_meta_d_mle` and
    :func:`metasignal.stdpy.trials_to_counts`.

    Parameters
    ----------
    d :
        Type-1 d-prime (sensitivity).
    noise :
        Standard deviation of metacognitive noise added to the type-1
        internal response for the type-2 judgment.
        Pass a two-element list ``[sigma_rS1, sigma_rS2]`` to simulate
        response-conditional noise.
    c :
        Type-1 decision criterion.
    nRatings :
        Number of confidence rating levels.
    nTrials :
        Number of trials to simulate.
    rng :
        Optional :class:`numpy.random.Generator` for reproducibility.
        If ``None``, uses :func:`numpy.random.default_rng`.

    Returns
    -------
    nR_S1, nR_S2 : ndarray
        Response-count vectors of length ``2 * nRatings``.

        *Ordering (same convention as metadpy / Maniscalco & Lau):*

        ``nR_S1[0 : nRatings]``   — S1 stimulus, S1 response, ratings high→low
        ``nR_S1[nRatings : 2*nRatings]`` — S1 stimulus, S2 response, ratings low→high

        Same layout for nR_S2.

    Examples
    --------
    >>> nR_S1, nR_S2 = type2_SDT_simuation(d=1.5, noise=0.0, c=0, nRatings=4, nTrials=1000)
    >>> len(nR_S1)
    8
    """
    if rng is None:
        rng = np.random.default_rng()

    # --- confidence criteria (same placement as metadpy) ---
    c1 = c + np.linspace(-1.5, -0.5, nRatings - 1)  # S1-response side (left of decision)
    c2 = c + np.linspace(0.5, 1.5, nRatings - 1)    # S2-response side (right of decision)

    # --- response-conditional noise ---
    if isinstance(noise, (list, tuple, np.ndarray)):
        response_conditional = True
        sigma_rS1 = float(noise[0])
        sigma_rS2 = float(noise[1])
    else:
        response_conditional = False
        sigma = float(noise)

    S1mu = -d / 2
    S2mu = d / 2

    # Counts: [nRatings] bins per response side
    nC_rS1 = np.zeros(nRatings)   # correct S1 responses (stim=S1, resp=S1)
    nI_rS1 = np.zeros(nRatings)   # incorrect S1 responses (stim=S2, resp=S1)
    nC_rS2 = np.zeros(nRatings)   # correct S2 responses (stim=S2, resp=S2)
    nI_rS2 = np.zeros(nRatings)   # incorrect S2 responses (stim=S1, resp=S2)

    for _ in range(nTrials):
        # Stimulus: 0 = S1, 1 = S2
        s = int(rng.integers(0, 2))

        # Type-1 internal response
        x = rng.normal(S2mu if s == 1 else S1mu, 1.0)

        # Metacognitive (type-2) response: add noise to x
        if response_conditional:
            sigma_now = sigma_rS1 if x < c else sigma_rS2
        else:
            sigma_now = sigma

        x2 = rng.normal(x, sigma_now) if sigma_now > 0 else x

        # Type-1 decision
        resp = int(x >= c)

        # Assign to confidence bin
        if s == 0 and resp == 0:          # S1 stim, S1 response (correct)
            i = np.where(np.append(c1, c) >= x2)[0]
            if len(i):
                nC_rS1[i.min()] += 1

        elif s == 0 and resp == 1:         # S1 stim, S2 response (incorrect)
            i = np.where(np.append(c, c2) <= x2)[0]
            if len(i):
                nI_rS2[i.max()] += 1

        elif s == 1 and resp == 0:         # S2 stim, S1 response (incorrect)
            i = np.where(np.append(c1, c) >= x2)[0]
            if len(i):
                nI_rS1[i.min()] += 1

        else:                              # S2 stim, S2 response (correct)
            i = np.where(np.append(c, c2) <= x2)[0]
            if len(i):
                nC_rS2[i.max()] += 1

    nR_S1 = np.concatenate([nC_rS1, nI_rS2])
    nR_S2 = np.concatenate([nI_rS1, nC_rS2])

    return nR_S1, nR_S2

type2_SDT_simuation_bayes(d=1.0, c=0.0, nRatings=4, nTrials=500, meta_noise=0.0, rng=None)

Simulate Type-2 SDT data using a Bayesian ideal observer.

Unlike :func:type2_SDT_simuation, which uses a Gaussian noise model (noisy access to the type-1 evidence signal), this function models an observer who bases confidence on the posterior probability of being correct — P(correct | x) — derived via Bayes' theorem.

A noiseless observer (meta_noise=0) is metacognitively ideal: meta-d' = d', giving M_ratio = 1 exactly. Adding meta_noise > 0 corrupts access to the posterior on the log-odds scale, reducing M_ratio below 1 in proportion to the noise.

This implements the framework described in Fleming & Dolan (2012) and separates the metacognitive noise model from the evidence noise model.

Parameters

d : Type-1 d-prime. c : Type-1 decision criterion. nRatings : Number of confidence rating levels. nTrials : Number of trials. meta_noise : Standard deviation of Gaussian noise added to the log-odds of P(correct|x) before binning. Set to 0 for an ideal Bayesian observer (M_ratio ≈ 1). Larger values reduce metacognitive efficiency. rng : Optional :class:numpy.random.Generator for reproducibility.

Returns

nR_S1, nR_S2 : ndarray Response-count vectors of length 2 * nRatings, in the same Maniscalco & Lau convention as :func:type2_SDT_simuation.

Notes

How confidence is computed

For evidence x drawn from N(±d/2, 1):

.. math::

P(S2|x) = \sigma(d \cdot x)   \quad (\text{equal priors, } c=0)

More generally with non-zero c we use the exact Gaussian ratio:

.. math::

P(S2|x) = \frac{\phi(x - d/2)}{\phi(x - d/2) + \phi(x + d/2)}

where :math:\phi is the standard normal PDF.

Confidence (posterior probability of being correct):

.. math::

p_{\text{correct}} = P(S2|x) \text{ if resp=S2, else } 1 - P(S2|x)

With meta-noise, noise is applied on the log-odds scale before converting back to probability:

.. math::

\ell_{\text{noisy}} = \text{logit}(p_{\text{correct}}) + \epsilon,
\quad \epsilon \sim N(0, \sigma_{\text{meta}})

Confidence binning

P(correct) ∈ [0.5, 1.0] is divided into nRatings equal-width bins. The bin edges are placed at 0.5 + k * 0.5/nRatings for k = 0, ..., nRatings. Bin index 0 = lowest confidence (P(correct) closest to 0.5), bin index nRatings-1 = highest confidence.

Examples

Ideal observer — expect M_ratio ≈ 1:

nR_S1, nR_S2 = type2_SDT_simuation_bayes(d=1.5, meta_noise=0.0, nTrials=5000)

Noisy observer — M_ratio drops below 1:

nR_S1, nR_S2 = type2_SDT_simuation_bayes(d=1.5, meta_noise=1.0, nTrials=5000)

Source code in src/metasignal/stdpy/simulate.py 
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
def type2_SDT_simuation_bayes(
    d: float = 1.0,
    c: float = 0.0,
    nRatings: int = 4,
    nTrials: int = 500,
    meta_noise: float = 0.0,
    rng: np.random.Generator | None = None,
) -> tuple[np.ndarray, np.ndarray]:
    """Simulate Type-2 SDT data using a Bayesian ideal observer.

    Unlike :func:`type2_SDT_simuation`, which uses a Gaussian noise model
    (noisy access to the type-1 evidence signal), this function models an
    observer who bases confidence on the **posterior probability of being
    correct** — P(correct | x) — derived via Bayes' theorem.

    A noiseless observer (``meta_noise=0``) is metacognitively ideal:
    meta-d' = d', giving M_ratio = 1 exactly.  Adding ``meta_noise > 0``
    corrupts access to the posterior on the log-odds scale, reducing M_ratio
    below 1 in proportion to the noise.

    This implements the framework described in Fleming & Dolan (2012) and
    separates the metacognitive noise model from the evidence noise model.

    Parameters
    ----------
    d :
        Type-1 d-prime.
    c :
        Type-1 decision criterion.
    nRatings :
        Number of confidence rating levels.
    nTrials :
        Number of trials.
    meta_noise :
        Standard deviation of Gaussian noise added to the **log-odds** of
        P(correct|x) before binning.  Set to 0 for an ideal Bayesian
        observer (M_ratio ≈ 1).  Larger values reduce metacognitive
        efficiency.
    rng :
        Optional :class:`numpy.random.Generator` for reproducibility.

    Returns
    -------
    nR_S1, nR_S2 : ndarray
        Response-count vectors of length ``2 * nRatings``, in the same
        Maniscalco & Lau convention as :func:`type2_SDT_simuation`.

    Notes
    -----
    **How confidence is computed**

    For evidence *x* drawn from N(±d/2, 1):

    .. math::

        P(S2|x) = \\sigma(d \\cdot x)   \\quad (\\text{equal priors, } c=0)

    More generally with non-zero *c* we use the exact Gaussian ratio:

    .. math::

        P(S2|x) = \\frac{\\phi(x - d/2)}{\\phi(x - d/2) + \\phi(x + d/2)}

    where :math:`\\phi` is the standard normal PDF.

    Confidence (posterior probability of being correct):

    .. math::

        p_{\\text{correct}} = P(S2|x) \\text{ if resp=S2, else } 1 - P(S2|x)

    With meta-noise, noise is applied on the log-odds scale before converting
    back to probability:

    .. math::

        \\ell_{\\text{noisy}} = \\text{logit}(p_{\\text{correct}}) + \\epsilon,
        \\quad \\epsilon \\sim N(0, \\sigma_{\\text{meta}})

    **Confidence binning**

    P(correct) ∈ [0.5, 1.0] is divided into *nRatings* equal-width bins.
    The bin edges are placed at ``0.5 + k * 0.5/nRatings`` for
    k = 0, ..., nRatings.  Bin index 0 = lowest confidence (P(correct)
    closest to 0.5), bin index nRatings-1 = highest confidence.

    Examples
    --------
    Ideal observer — expect M_ratio ≈ 1:

    >>> nR_S1, nR_S2 = type2_SDT_simuation_bayes(d=1.5, meta_noise=0.0, nTrials=5000)

    Noisy observer — M_ratio drops below 1:

    >>> nR_S1, nR_S2 = type2_SDT_simuation_bayes(d=1.5, meta_noise=1.0, nTrials=5000)
    """
    if rng is None:
        rng = np.random.default_rng()

    S1mu = -d / 2
    S2mu = d / 2

    nC_rS1 = np.zeros(nRatings)
    nI_rS1 = np.zeros(nRatings)
    nC_rS2 = np.zeros(nRatings)
    nI_rS2 = np.zeros(nRatings)

    # Equal-width bin edges on [0.5, 1.0]
    bin_edges = np.linspace(0.5, 1.0, nRatings + 1)

    for _ in range(nTrials):
        s = int(rng.integers(0, 2))
        x = rng.normal(S2mu if s == 1 else S1mu, 1.0)
        resp = int(x >= c)

        # --- Bayesian posterior P(S2|x) via Gaussian ratio ---
        log_lik_s2 = norm.logpdf(x, S2mu, 1.0)
        log_lik_s1 = norm.logpdf(x, S1mu, 1.0)
        # log P(S2|x) - log P(S1|x) = log_lik_s2 - log_lik_s1  (equal priors)
        log_odds_s2 = log_lik_s2 - log_lik_s1
        p_s2 = float(expit(log_odds_s2))

        # P(correct | x, response)
        p_correct = p_s2 if resp == 1 else 1.0 - p_s2

        # --- optional log-odds noise on the metacognitive posterior ---
        if meta_noise > 0:
            # logit of p_correct, clipped to avoid ±inf
            p_c_clipped = float(np.clip(p_correct, 1e-9, 1 - 1e-9))
            log_odds_c = np.log(p_c_clipped / (1.0 - p_c_clipped))
            log_odds_noisy = log_odds_c + rng.normal(0.0, meta_noise)
            p_correct = float(expit(log_odds_noisy))

        # --- bin P(correct) into 0..nRatings-1 ---
        # np.digitize returns 1-based; clip to [0, nRatings-1]
        conf_bin = int(np.clip(np.digitize(p_correct, bin_edges) - 1, 0, nRatings - 1))

        # --- assign to count arrays ---
        # S1-response counts: index 0 = most confident (mirrors evidence model convention)
        # S2-response counts: index nRatings-1 = most confident
        if s == 0 and resp == 0:        # S1 stim, S1 resp (correct) — most conf = bin 0
            nC_rS1[nRatings - 1 - conf_bin] += 1
        elif s == 0 and resp == 1:      # S1 stim, S2 resp (incorrect)
            nI_rS2[conf_bin] += 1
        elif s == 1 and resp == 0:      # S2 stim, S1 resp (incorrect) — most conf = bin 0
            nI_rS1[nRatings - 1 - conf_bin] += 1
        else:                           # S2 stim, S2 resp (correct)
            nC_rS2[conf_bin] += 1

    nR_S1 = np.concatenate([nC_rS1, nI_rS2])
    nR_S2 = np.concatenate([nI_rS1, nC_rS2])

    return nR_S1, nR_S2

ratings2df(nR_S1, nR_S2)

Convert nR_S1 / nR_S2 count arrays to a trial-level DataFrame.

Port of metadpy.utils.ratings2df.

Parameters

nR_S1, nR_S2 : Count vectors, length 2 * nRatings.

Returns

df : pd.DataFrame Columns: Stimuli (0/1), Responses (0/1), Accuracy (0/1), Confidence (1..nRatings), nTrial. Rows are shuffled.

Source code in src/metasignal/stdpy/simulate.py 
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
def ratings2df(nR_S1: np.ndarray, nR_S2: np.ndarray) -> "pd.DataFrame":
    """Convert nR_S1 / nR_S2 count arrays to a trial-level DataFrame.

    Port of ``metadpy.utils.ratings2df``.

    Parameters
    ----------
    nR_S1, nR_S2 :
        Count vectors, length ``2 * nRatings``.

    Returns
    -------
    df : pd.DataFrame
        Columns: ``Stimuli`` (0/1), ``Responses`` (0/1), ``Accuracy`` (0/1),
        ``Confidence`` (1..nRatings), ``nTrial``. Rows are shuffled.
    """
    import pandas as _pd
    nR_S1 = np.asarray(nR_S1, dtype=int)
    nR_S2 = np.asarray(nR_S2, dtype=int)
    nRatings = len(nR_S1) // 2

    rows: list = []
    for i in range(nRatings):
        conf = nRatings - i           # high confidence at low index
        for _ in range(nR_S1[i]):
            rows.append({"Stimuli": 0, "Responses": 0, "Accuracy": 1, "Confidence": conf})
        for _ in range(nR_S2[i]):
            rows.append({"Stimuli": 1, "Responses": 0, "Accuracy": 0, "Confidence": conf})
        for _ in range(nR_S1[nRatings + i]):
            rows.append({"Stimuli": 0, "Responses": 1, "Accuracy": 0, "Confidence": i + 1})
        for _ in range(nR_S2[nRatings + i]):
            rows.append({"Stimuli": 1, "Responses": 1, "Accuracy": 1, "Confidence": i + 1})

    df = _pd.DataFrame(rows).sample(frac=1).reset_index(drop=True)
    df["nTrial"] = np.arange(len(df))
    return df

trialSimulation(d=1.0, metad=2.0, c=0.0, nRatings=4, nTrials=500, rng=None)

Simulate a single participant using closed-form SDT probabilities.

Computes exact multinomial probabilities from normal-CDF integrals and draws counts — much faster than trial-by-trial noise simulation. Port of metadpy.utils.trialSimulation.

Parameters

d : Type-1 d-prime. metad : Type-2 meta-d' in d-prime units. c : Type-1 decision criterion. nRatings : Number of confidence rating levels. nTrials : Number of trials. rng : Optional :class:numpy.random.Generator for reproducibility.

Returns

df : pd.DataFrame Columns: Stimuli, Responses, Accuracy, Confidence, nTrial.

Source code in src/metasignal/stdpy/simulate.py 
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
def trialSimulation(
    d: float = 1.0,
    metad: float = 2.0,
    c: float = 0.0,
    nRatings: int = 4,
    nTrials: int = 500,
    rng: "np.random.Generator | None" = None,
) -> "pd.DataFrame":
    """Simulate a single participant using closed-form SDT probabilities.

    Computes exact multinomial probabilities from normal-CDF integrals and
    draws counts — much faster than trial-by-trial noise simulation.
    Port of ``metadpy.utils.trialSimulation``.

    Parameters
    ----------
    d :
        Type-1 d-prime.
    metad :
        Type-2 meta-d' in d-prime units.
    c :
        Type-1 decision criterion.
    nRatings :
        Number of confidence rating levels.
    nTrials :
        Number of trials.
    rng :
        Optional :class:`numpy.random.Generator` for reproducibility.

    Returns
    -------
    df : pd.DataFrame
        Columns: ``Stimuli``, ``Responses``, ``Accuracy``, ``Confidence``, ``nTrial``.
    """
    if rng is None:
        rng = np.random.default_rng()

    c1 = c + np.linspace(-1.5, -0.5, nRatings - 1)
    c2 = c + np.linspace(0.5,  1.5,  nRatings - 1)

    H  = round((1 - norm.cdf(c,  d / 2)) * (nTrials / 2))
    FA = round((1 - norm.cdf(c, -d / 2)) * (nTrials / 2))
    CR = round(norm.cdf(c, -d / 2)       * (nTrials / 2))
    M  = round(norm.cdf(c,  d / 2)       * (nTrials / 2))

    S1mu = -metad / 2
    S2mu =  metad / 2

    C_area_rS1 = norm.cdf(c, S1mu)
    I_area_rS1 = norm.cdf(c, S2mu)
    C_area_rS2 = 1 - norm.cdf(c, S2mu)
    I_area_rS2 = 1 - norm.cdf(c, S1mu)

    t2c1x = np.concatenate([[-np.inf], c1, [c], c2, [np.inf]])

    prC_rS1, prI_rS1, prC_rS2, prI_rS2 = [], [], [], []
    for i in range(nRatings):
        prC_rS1.append((norm.cdf(t2c1x[i + 1], S1mu) - norm.cdf(t2c1x[i], S1mu)) / C_area_rS1)
        prI_rS1.append((norm.cdf(t2c1x[i + 1], S2mu) - norm.cdf(t2c1x[i], S2mu)) / I_area_rS1)
        prC_rS2.append(
            ((1 - norm.cdf(t2c1x[nRatings + i],     S2mu))
           - (1 - norm.cdf(t2c1x[nRatings + i + 1], S2mu))) / C_area_rS2
        )
        prI_rS2.append(
            ((1 - norm.cdf(t2c1x[nRatings + i],     S1mu))
           - (1 - norm.cdf(t2c1x[nRatings + i + 1], S1mu))) / I_area_rS2
        )

    def _norm(p: list) -> np.ndarray:
        a = np.clip(np.array(p, dtype=float), 0, None)
        return a / a.sum()

    nC_rS1 = rng.multinomial(CR, _norm(prC_rS1))
    nI_rS1 = rng.multinomial(M,  _norm(prI_rS1))
    nC_rS2 = rng.multinomial(H,  _norm(prC_rS2))
    nI_rS2 = rng.multinomial(FA, _norm(prI_rS2))

    nr_s1 = np.concatenate([nC_rS1, nI_rS2])
    nr_s2 = np.concatenate([nI_rS1, nC_rS2])

    return ratings2df(nr_s1, nr_s2)

responseSimulation(d=1.0, metad=2.0, c=0.0, nRatings=4, nTrials=500, nSubjects=1, rng=None)

Simulate responses for one or more participants.

Calls :func:trialSimulation once per subject and stacks the results. Port of metadpy.utils.responseSimulation.

Parameters

d : Type-1 d-prime (same for all subjects). metad : Meta-d' (same for all subjects). c : Type-1 criterion (same for all subjects). nRatings : Number of confidence rating levels. nTrials : Trials per subject. nSubjects : Number of subjects. rng : Optional :class:numpy.random.Generator.

Returns

df : pd.DataFrame Columns: Stimuli, Responses, Accuracy, Confidence, nTrial, Subject.

Source code in src/metasignal/stdpy/simulate.py 
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
def responseSimulation(
    d: float = 1.0,
    metad: float = 2.0,
    c: float = 0.0,
    nRatings: int = 4,
    nTrials: int = 500,
    nSubjects: int = 1,
    rng: "np.random.Generator | None" = None,
) -> "pd.DataFrame":
    """Simulate responses for one or more participants.

    Calls :func:`trialSimulation` once per subject and stacks the results.
    Port of ``metadpy.utils.responseSimulation``.

    Parameters
    ----------
    d :
        Type-1 d-prime (same for all subjects).
    metad :
        Meta-d' (same for all subjects).
    c :
        Type-1 criterion (same for all subjects).
    nRatings :
        Number of confidence rating levels.
    nTrials :
        Trials per subject.
    nSubjects :
        Number of subjects.
    rng :
        Optional :class:`numpy.random.Generator`.

    Returns
    -------
    df : pd.DataFrame
        Columns: ``Stimuli``, ``Responses``, ``Accuracy``, ``Confidence``,
        ``nTrial``, ``Subject``.
    """
    import pandas as _pd
    if rng is None:
        rng = np.random.default_rng()

    frames = []
    for sub in range(nSubjects):
        df = trialSimulation(d=d, metad=metad, c=c,
                             nRatings=nRatings, nTrials=nTrials, rng=rng)
        df["Subject"] = sub
        frames.append(df)

    return _pd.concat(frames, ignore_index=True)

pairedResponseSimulation(d=1.0, d_sigma=0.1, mRatio=None, mRatio_sigma=0.2, mRatio_rho=0.0, c=0.0, c_sigma=0.1, nRatings=4, nTrials=500, nSubjects=20, rng=None)

Simulate a two-condition repeated-measures group design.

Draws correlated M-ratio pairs per subject from a bivariate normal distribution, then simulates both conditions via :func:trialSimulation. Port of metadpy.utils.pairedResponseSimulation.

Parameters

d : Group-mean type-1 d-prime. d_sigma : Between-subject SD for d-prime. mRatio : [mRatio_cond1, mRatio_cond2]. Defaults to [1.0, 0.6]. mRatio_sigma : SD of M-ratio across subjects (applied to both conditions). mRatio_rho : Correlation between the two per-subject M-ratios. c : Group-mean type-1 criterion. c_sigma : Between-subject SD for criterion. nRatings : Number of confidence rating levels. nTrials : Trials per subject per condition. nSubjects : Number of subjects. rng : Optional :class:numpy.random.Generator.

Returns

df : pd.DataFrame Columns: Stimuli, Responses, Accuracy, Confidence, nTrial, Subject, Condition.

Source code in src/metasignal/stdpy/simulate.py 
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
def pairedResponseSimulation(
    d: float = 1.0,
    d_sigma: float = 0.1,
    mRatio=None,
    mRatio_sigma: float = 0.2,
    mRatio_rho: float = 0.0,
    c: float = 0.0,
    c_sigma: float = 0.1,
    nRatings: int = 4,
    nTrials: int = 500,
    nSubjects: int = 20,
    rng: "np.random.Generator | None" = None,
) -> "pd.DataFrame":
    """Simulate a two-condition repeated-measures group design.

    Draws correlated M-ratio pairs per subject from a bivariate normal
    distribution, then simulates both conditions via :func:`trialSimulation`.
    Port of ``metadpy.utils.pairedResponseSimulation``.

    Parameters
    ----------
    d :
        Group-mean type-1 d-prime.
    d_sigma :
        Between-subject SD for d-prime.
    mRatio :
        ``[mRatio_cond1, mRatio_cond2]``. Defaults to ``[1.0, 0.6]``.
    mRatio_sigma :
        SD of M-ratio across subjects (applied to both conditions).
    mRatio_rho :
        Correlation between the two per-subject M-ratios.
    c :
        Group-mean type-1 criterion.
    c_sigma :
        Between-subject SD for criterion.
    nRatings :
        Number of confidence rating levels.
    nTrials :
        Trials per subject per condition.
    nSubjects :
        Number of subjects.
    rng :
        Optional :class:`numpy.random.Generator`.

    Returns
    -------
    df : pd.DataFrame
        Columns: ``Stimuli``, ``Responses``, ``Accuracy``, ``Confidence``,
        ``nTrial``, ``Subject``, ``Condition``.
    """
    import pandas as _pd
    if mRatio is None:
        mRatio = [1.0, 0.6]
    if rng is None:
        rng = np.random.default_rng()

    cov = np.array([
        [mRatio_sigma ** 2, mRatio_rho * mRatio_sigma ** 2],
        [mRatio_rho * mRatio_sigma ** 2, mRatio_sigma ** 2],
    ])

    frames = []
    for sub in range(nSubjects):
        mr_pair = rng.multivariate_normal(mRatio, cov)
        for cond in range(2):
            d_sub     = float(rng.normal(d, d_sigma))
            c_sub     = float(rng.normal(c, c_sigma))
            metad_sub = float(mr_pair[cond] * d_sub)
            df = trialSimulation(d=d_sub, metad=metad_sub, c=c_sub,
                                 nRatings=nRatings, nTrials=nTrials, rng=rng)
            df["Subject"]   = sub
            df["Condition"] = cond
            frames.append(df)

    return _pd.concat(frames, ignore_index=True)

discreteRatings(ratings, nbins=4, verbose=True, ignore_invalid=False)

Bin continuous confidence ratings into discrete levels via quantiles.

Handles floor/ceiling biases by treating extreme values as their own bin after re-estimating inner quantile boundaries from the remaining data. Port of metadpy.utils.discreteRatings.

Parameters

ratings : Continuous confidence ratings (any numeric scale). nbins : Target number of discrete bins. verbose : Print a message when bias correction is applied. ignore_invalid : If False, raise ValueError when both floor and ceiling are biased. Set True to proceed anyway.

Returns

discrete : np.ndarray Integer ratings from 1 to nbins. info : dict 'confBins', 'rebin' ([1] if corrected, else [0]), 'binCount' (per-bin trial counts).

Source code in src/metasignal/stdpy/simulate.py 
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
def discreteRatings(ratings, nbins: int = 4, verbose: bool = True,
                    ignore_invalid: bool = False):
    """Bin continuous confidence ratings into discrete levels via quantiles.

    Handles floor/ceiling biases by treating extreme values as their own bin
    after re-estimating inner quantile boundaries from the remaining data.
    Port of ``metadpy.utils.discreteRatings``.

    Parameters
    ----------
    ratings :
        Continuous confidence ratings (any numeric scale).
    nbins :
        Target number of discrete bins.
    verbose :
        Print a message when bias correction is applied.
    ignore_invalid :
        If ``False``, raise ``ValueError`` when both floor and ceiling are
        biased. Set ``True`` to proceed anyway.

    Returns
    -------
    discrete : np.ndarray
        Integer ratings from 1 to ``nbins``.
    info : dict
        ``'confBins'``, ``'rebin'`` (``[1]`` if corrected, else ``[0]``),
        ``'binCount'`` (per-bin trial counts).
    """
    ratings = np.asarray(ratings, dtype=float)
    info: dict = {}
    masks: list = []

    conf_bins    = np.quantile(ratings, np.linspace(0, 1, nbins + 1))
    floor_bias   = conf_bins[0] == conf_bins[1]
    ceiling_bias = conf_bins[nbins - 1] == conf_bins[nbins]

    if floor_bias and ceiling_bias:
        if not ignore_invalid:
            raise ValueError(
                "Rating scale has too many identical extreme values to discretise. "
                "Pass ignore_invalid=True to proceed."
            )
        if verbose:
            print("Correcting for bias in both low and high confidence ratings.")
        lo = conf_bins[0]
        hi = conf_bins[-1]
        inner_ratings = ratings[(ratings != lo) & (ratings != hi)]
        n_inner_bins = max(nbins - 2, 0)
        masks.append(ratings == lo)
        if n_inner_bins > 0 and len(inner_ratings) > 0:
            inner = np.quantile(inner_ratings, np.linspace(0, 1, n_inner_bins + 1))
            for b in range(n_inner_bins):
                masks.append(
                    (ratings >= inner[b]) & (ratings <= inner[b + 1])
                    & (ratings != lo) & (ratings != hi)
                )
        masks.append(ratings == hi)
        info["confBins"] = [lo, hi]
        info["rebin"]    = [1]

    elif ceiling_bias:
        if verbose:
            print("Correcting for bias in high confidence ratings.")
        hi = conf_bins[-1]
        inner = np.quantile(ratings[ratings != hi], np.linspace(0, 1, nbins))
        for b in range(len(inner) - 1):
            masks.append((ratings >= inner[b]) & (ratings <= inner[b + 1]))
        masks.append(ratings == hi)
        info["confBins"] = [inner, hi]
        info["rebin"]    = [1]

    elif floor_bias:
        if verbose:
            print("Correcting for bias in low confidence ratings.")
        lo = conf_bins[1]
        masks.append(ratings == lo)
        inner = np.quantile(ratings[ratings != lo], np.linspace(0, 1, nbins))
        for b in range(1, len(inner)):
            masks.append((ratings >= inner[b - 1]) & (ratings <= inner[b]))
        info["confBins"] = [lo, inner]
        info["rebin"]    = [1]

    else:
        for b in range(len(conf_bins) - 1):
            masks.append((ratings >= conf_bins[b]) & (ratings <= conf_bins[b + 1]))
        info["confBins"] = conf_bins
        info["rebin"]    = [0]

    discrete = np.zeros(len(ratings), dtype=int)
    for b, mask in enumerate(masks):
        discrete[mask] = b
    discrete += 1  # 1-based

    info["binCount"] = [int(m.sum()) for m in masks]
    return discrete, info

analysis — Statistical inference

The analysis sub-package provides non-parametric inference tools that operate over the 26-element array (20 measures + 6 meta-d' model-fit diagnostics) produced by compute_all_measures.

bootstrap_measure(stim, resp, conf, n_ratings, measure_index, n_boot=2000, ci=0.95, rng=None)

Return a bootstrap confidence interval for one element of the 26-measure array.

Resamples trials with replacement n_boot times and computes the requested measure on each resample. The interval is percentile-based.

Parameters:

Name Type Description Default
stim ndarray

Stimulus array (binary, same encoding as compute_all_measures).

 required
resp ndarray

Response array (binary).

 required
conf ndarray

Confidence rating array (1 to n_ratings).

 required
n_ratings int

Number of confidence rating categories.

 required
measure_index int

Index into the 26-element output of compute_all_measures. See that function's docstring for the index→measure mapping.

 required
n_boot int

Number of bootstrap resamples. Default 2000.

 2000
ci float

Coverage of the interval, e.g. 0.95 for a 95% CI. Default 0.95.

 0.95
rng Generator | None

Optional numpy.random.Generator for reproducibility.

 None

Returns:

Type Description
tuple[float, float]

Tuple (lower, upper) confidence bounds.

Raises:

Type Description
ValueError

If measure_index is outside [0, 25].
Source code in src/metasignal/analysis/bootstrap.py 
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
def bootstrap_measure(
    stim: np.ndarray,
    resp: np.ndarray,
    conf: np.ndarray,
    n_ratings: int,
    measure_index: int,
    n_boot: int = 2000,
    ci: float = 0.95,
    rng: np.random.Generator | None = None,
) -> tuple[float, float]:
    """Return a bootstrap confidence interval for one element of the 26-measure array.

    Resamples trials with replacement ``n_boot`` times and computes the
    requested measure on each resample. The interval is percentile-based.

    Args:
        stim: Stimulus array (binary, same encoding as ``compute_all_measures``).
        resp: Response array (binary).
        conf: Confidence rating array (1 to n_ratings).
        n_ratings: Number of confidence rating categories.
        measure_index: Index into the 26-element output of ``compute_all_measures``.
            See that function's docstring for the index→measure mapping.
        n_boot: Number of bootstrap resamples. Default 2000.
        ci: Coverage of the interval, e.g. 0.95 for a 95% CI. Default 0.95.
        rng: Optional ``numpy.random.Generator`` for reproducibility.

    Returns:
        Tuple ``(lower, upper)`` confidence bounds.

    Raises:
        ValueError: If ``measure_index`` is outside [0, 25].
    """
    if not 0 <= measure_index <= 25:
        msg = f"measure_index must be in [0, 25], got {measure_index}"
        raise ValueError(msg)

    if rng is None:
        rng = np.random.default_rng()

    stim = np.asarray(stim)
    resp = np.asarray(resp)
    conf = np.asarray(conf)
    n = len(stim)

    boot_vals: list[float] = []
    for _ in range(n_boot):
        idx = rng.integers(0, n, size=n)
        result = compute_all_measures(stim[idx], resp[idx], conf[idx], n_ratings)
        val = result[measure_index]
        if not np.isnan(val):
            boot_vals.append(float(val))

    if len(boot_vals) == 0:
        warnings.warn(
            f"bootstrap_measure: all {n_boot} resamples produced NaN for "
            f"measure_index={measure_index}; returning (nan, nan).",
            stacklevel=2,
        )
        return (float("nan"), float("nan"))

    if len(boot_vals) < 0.5 * n_boot:
        warnings.warn(
            f"bootstrap_measure: only {len(boot_vals)} of {n_boot} requested "
            f"resamples produced a finite value for measure_index={measure_index} "
            "(the rest were NaN, e.g. due to MLE non-convergence); the returned "
            "CI is based on this smaller effective sample.",
            stacklevel=2,
        )

    alpha = 1.0 - ci
    lo = float(np.percentile(boot_vals, 100 * alpha / 2))
    hi = float(np.percentile(boot_vals, 100 * (1 - alpha / 2)))
    return lo, hi

permutation_test(stim_a, resp_a, conf_a, stim_b, resp_b, conf_b, n_ratings, measure_index, n_perm=5000, rng=None)

Two-sided permutation test for a difference in one SDT measure between conditions.

Shuffles trial condition labels n_perm times and computes the difference in the measure on each shuffle to build a null distribution. Returns the two-sided p-value and the observed difference.

Parameters:

Name Type Description Default
stim_a ndarray

Stimulus array for condition A.

 required
resp_a ndarray

Response array for condition A.

 required
conf_a ndarray

Confidence array for condition A.

 required
stim_b ndarray

Stimulus array for condition B.

 required
resp_b ndarray

Response array for condition B.

 required
conf_b ndarray

Confidence array for condition B.

 required
n_ratings int

Number of confidence rating categories (shared across conditions).

 required
measure_index int

Index into the 26-element output of compute_all_measures.

 required
n_perm int

Number of permutations. Default 5000.

 5000
rng Generator | None

Optional numpy.random.Generator for reproducibility.

 None

Returns:

Type Description
float

Tuple (p_value, observed_difference). observed_difference is
float

measure(A) - measure(B).

Raises:

Type Description
ValueError

If measure_index is outside [0, 25].
Source code in src/metasignal/analysis/permutation.py 
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
def permutation_test(
    stim_a: np.ndarray,
    resp_a: np.ndarray,
    conf_a: np.ndarray,
    stim_b: np.ndarray,
    resp_b: np.ndarray,
    conf_b: np.ndarray,
    n_ratings: int,
    measure_index: int,
    n_perm: int = 5000,
    rng: np.random.Generator | None = None,
) -> tuple[float, float]:
    """Two-sided permutation test for a difference in one SDT measure between conditions.

    Shuffles trial condition labels ``n_perm`` times and computes the
    difference in the measure on each shuffle to build a null distribution.
    Returns the two-sided p-value and the observed difference.

    Args:
        stim_a: Stimulus array for condition A.
        resp_a: Response array for condition A.
        conf_a: Confidence array for condition A.
        stim_b: Stimulus array for condition B.
        resp_b: Response array for condition B.
        conf_b: Confidence array for condition B.
        n_ratings: Number of confidence rating categories (shared across conditions).
        measure_index: Index into the 26-element output of ``compute_all_measures``.
        n_perm: Number of permutations. Default 5000.
        rng: Optional ``numpy.random.Generator`` for reproducibility.

    Returns:
        Tuple ``(p_value, observed_difference)``.  ``observed_difference`` is
        ``measure(A) - measure(B)``.

    Raises:
        ValueError: If ``measure_index`` is outside [0, 25].
    """
    if not 0 <= measure_index <= 25:
        msg = f"measure_index must be in [0, 25], got {measure_index}"
        raise ValueError(msg)

    if rng is None:
        rng = np.random.default_rng()

    stim_a, resp_a, conf_a = (np.asarray(x) for x in (stim_a, resp_a, conf_a))
    stim_b, resp_b, conf_b = (np.asarray(x) for x in (stim_b, resp_b, conf_b))

    def _measure(s: np.ndarray, r: np.ndarray, c: np.ndarray) -> float:
        return float(compute_all_measures(s, r, c, n_ratings)[measure_index])

    obs_a = _measure(stim_a, resp_a, conf_a)
    obs_b = _measure(stim_b, resp_b, conf_b)
    obs_diff = obs_a - obs_b

    # Pool all trials for shuffling
    stim_all = np.concatenate([stim_a, stim_b])
    resp_all = np.concatenate([resp_a, resp_b])
    conf_all = np.concatenate([conf_a, conf_b])
    n_a = len(stim_a)
    n_total = len(stim_all)

    null_diffs: list[float] = []
    for _ in range(n_perm):
        perm = rng.permutation(n_total)
        idx_a, idx_b = perm[:n_a], perm[n_a:]
        val_a = _measure(stim_all[idx_a], resp_all[idx_a], conf_all[idx_a])
        val_b = _measure(stim_all[idx_b], resp_all[idx_b], conf_all[idx_b])
        diff = val_a - val_b
        if not np.isnan(diff):
            null_diffs.append(diff)

    if len(null_diffs) == 0:
        warnings.warn(
            f"permutation_test: all {n_perm} permutations produced NaN for "
            f"measure_index={measure_index}; returning (nan, observed_difference).",
            stacklevel=2,
        )
        return (float("nan"), obs_diff)

    if len(null_diffs) < 0.5 * n_perm:
        warnings.warn(
            f"permutation_test: only {len(null_diffs)} of {n_perm} requested "
            f"permutations produced a finite value for measure_index={measure_index} "
            "(the rest were NaN); the returned p-value is based on this smaller "
            "effective null distribution.",
            stacklevel=2,
        )

    null = np.array(null_diffs)
    n_extreme = int(np.sum(np.abs(null) >= abs(obs_diff)))
    p_value = (n_extreme + 1) / (len(null) + 1)
    return p_value, obs_diff

group_summary(participants, n_ratings)

Compute group-level descriptive statistics over the 26-measure array.

Runs compute_all_measures for each participant and aggregates across the group, ignoring NaN values (participants where a measure could not be estimated are excluded from that measure's statistics).

Parameters:

Name Type Description Default
participants list[tuple[ndarray, ndarray, ndarray]]

List of (stim, resp, conf) tuples, one per participant.

 required
n_ratings int

Number of confidence rating categories.

 required

Returns:

Type Description
dict[str, Any]

Dictionary with keys:
dict[str, Any]
  • "individual" — np.ndarray of shape (n_participants, 26)
dict[str, Any]
  • "mean" — np.ndarray of shape (26,), nanmean across participants
dict[str, Any]
  • "median" — np.ndarray of shape (26,)
dict[str, Any]
  • "sem" — np.ndarray of shape (26,), standard error of the mean
dict[str, Any]
  • "n_valid" — np.ndarray of shape (26,) int, participants with non-NaN values per measure
dict[str, Any]
  • "labels" — list of measure name strings (length 26)
Source code in src/metasignal/analysis/group.py 
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
def group_summary(
    participants: list[tuple[np.ndarray, np.ndarray, np.ndarray]],
    n_ratings: int,
) -> dict[str, Any]:
    """Compute group-level descriptive statistics over the 26-measure array.

    Runs ``compute_all_measures`` for each participant and aggregates across
    the group, ignoring NaN values (participants where a measure could not
    be estimated are excluded from that measure's statistics).

    Args:
        participants: List of ``(stim, resp, conf)`` tuples, one per participant.
        n_ratings: Number of confidence rating categories.

    Returns:
        Dictionary with keys:

        - ``"individual"`` — ``np.ndarray`` of shape ``(n_participants, 26)``
        - ``"mean"`` — ``np.ndarray`` of shape ``(26,)``, nanmean across participants
        - ``"median"`` — ``np.ndarray`` of shape ``(26,)``
        - ``"sem"`` — ``np.ndarray`` of shape ``(26,)``, standard error of the mean
        - ``"n_valid"`` — ``np.ndarray`` of shape ``(26,)`` int, participants with
            non-NaN values per measure
        - ``"labels"`` — list of measure name strings (length 26)
    """
    n_participants = len(participants)
    individual = np.full((n_participants, N_MEASURES), np.nan)

    for i, (stim, resp, conf) in enumerate(participants):
        individual[i] = compute_all_measures(
            np.asarray(stim), np.asarray(resp), np.asarray(conf), n_ratings
        )

    n_valid = np.sum(~np.isnan(individual), axis=0).astype(int)
    mean = np.nanmean(individual, axis=0)
    median = np.nanmedian(individual, axis=0)
    # Require at least 2 valid values for a meaningful SEM; return NaN otherwise
    std = np.nanstd(individual, axis=0, ddof=1)
    sem = np.where(n_valid >= 2, std / np.sqrt(n_valid), np.nan)

    return {
        "individual": individual,
        "mean": mean,
        "median": median,
        "sem": sem,
        "n_valid": n_valid,
        "labels": MEASURE_LABELS,
    }

itmc — Information-theoretic metacognition (experimental)

The itmc sub-package implements the information-theoretic metacognition framework of Dayan (2023) — Metacognitive Information Theory, Open Mind, 7, 392–411, doi:10.1162/opmi_a_00091 — measuring metacognitive sensitivity as mutual information between accuracy and confidence. meta_I, meta_Ir1, and meta_Ir2 are Dayan's own proposed measures; meta_Ir1_acc and RMI are related measures introduced by Rausch et al. (2025) — statConfR: An R Package for Static Models of Decision Confidence and Metacognition, Journal of Open Source Software, 10(106), 6966, doi:10.21105/joss.06966 — that extend Dayan's framework but aren't discussed in his original paper. It is a pre-1.0 component and its API may still change between releases. Its test suite (tests/test_itmc.py) covers bounds, monotonicity, and internal-consistency properties; backend='statconfr' is additionally cross-validated against the real R statConfR package it was ported from, and meta_I independently reproduces Dayan (2023)'s own hand-worked numerical example exactly (0.0943 vs. his published 0.094) — see analysis/itmc_comparison/ANALYSIS_REPORT.md for the full agreement report.

meta_I(stimulus, response, rating, *, backend='simple', bias_correction=False, seed=42)

Mutual information between confidence and accuracy (bits).

Parameters

stimulus, response, rating: Trial-level arrays. stimulus/response are 0/1; rating is 1…n_ratings. backend: 'simple': H₂(acc) − H₂(acc|conf). Fast, directly measures how much confidence predicts correctness.

``'statconfr'``: I(stimulus; graded_response) − I_min(prior, acc).
Builds the full 2×2K joint table, computes Shannon MI, then subtracts
the analytic minimum achievable at this accuracy.  Matches the R package
exactly.

bias_correction: Subtract the estimated positive sampling bias. 'simple': permutation of accuracy labels (2000 shuffles). 'statconfr': multinomial resampling of each table row (1000 sims), matching statConfR's bias_reduction=TRUE. Returns mean(simulated), which can be higher or lower than the observed value. seed: RNG seed for bias reduction.

Source code in src/metasignal/itmc/measures.py 
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
def meta_I(
    stimulus: np.ndarray,
    response: np.ndarray,
    rating: np.ndarray,
    *,
    backend: Backend = "simple",
    bias_correction: bool = False,
    seed: int = 42,
) -> float:
    """Mutual information between confidence and accuracy (bits).

    Parameters
    ----------
    stimulus, response, rating:
        Trial-level arrays. stimulus/response are 0/1; rating is 1…n_ratings.
    backend:
        ``'simple'``: H₂(acc) − H₂(acc|conf).  Fast, directly measures how
        much confidence predicts correctness.

        ``'statconfr'``: I(stimulus; graded_response) − I_min(prior, acc).
        Builds the full 2×2K joint table, computes Shannon MI, then subtracts
        the analytic minimum achievable at this accuracy.  Matches the R package
        exactly.
    bias_correction:
        Subtract the estimated positive sampling bias.
        ``'simple'``: permutation of accuracy labels (2000 shuffles).
        ``'statconfr'``: multinomial resampling of each table row (1000 sims),
        matching statConfR's ``bias_reduction=TRUE``.  Returns mean(simulated),
        which can be higher or lower than the observed value.
    seed:
        RNG seed for bias reduction.
    """
    stimulus = np.asarray(stimulus)
    response = np.asarray(response)
    rating   = np.asarray(rating)

    if backend == "simple":
        acc = (stimulus == response).astype(int)
        mi  = _simple_meta_I(acc, rating)
        if bias_correction:
            mi -= _simple_bias(acc, rating, seed=seed)
        return float(mi)

    if backend == "statconfr":
        table = _build_contingency_table(stimulus, response, rating)
        if bias_correction:
            bias = _statconfr_bias(table, _sc_meta_I, seed=seed)
            return float(_sc_meta_I(table) - bias)
        return _sc_meta_I(table)

    raise ValueError(f"Unknown backend {backend!r}. Choose 'simple' or 'statconfr'.")

meta_Ir1(stimulus, response, rating, dprime=None, *, backend='simple', bias_correction=False, seed=42)

Gaussian-normalised relative metacognitive efficiency (meta-I₁ʳ).

meta-I₁ʳ = meta-I / meta-I_Gaussian(d′)

Parameters

dprime: First-order sensitivity. If None, estimated from stimulus/response. backend: 'simple': Gaussian reference via Monte Carlo simulation with equal-count confidence bins.

``'statconfr'``: Gaussian reference via continuous 5000-point grid
integration of the ideal Gaussian classifier, matching the R package.
Source code in src/metasignal/itmc/measures.py 
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
def meta_Ir1(
    stimulus: np.ndarray,
    response: np.ndarray,
    rating: np.ndarray,
    dprime: float | None = None,
    *,
    backend: Backend = "simple",
    bias_correction: bool = False,
    seed: int = 42,
) -> float:
    """Gaussian-normalised relative metacognitive efficiency (meta-I₁ʳ).

    meta-I₁ʳ = meta-I / meta-I_Gaussian(d′)

    Parameters
    ----------
    dprime:
        First-order sensitivity. If None, estimated from stimulus/response.
    backend:
        ``'simple'``: Gaussian reference via Monte Carlo simulation with
        equal-count confidence bins.

        ``'statconfr'``: Gaussian reference via continuous 5000-point grid
        integration of the ideal Gaussian classifier, matching the R package.
    """
    stimulus = np.asarray(stimulus)
    response = np.asarray(response)
    rating   = np.asarray(rating)

    if backend == "simple":
        from metasignal.stdpy.core import compute_sdt_resp
        acc = (stimulus == response).astype(int)
        mi  = _simple_meta_I(acc, rating)
        if bias_correction:
            mi -= _simple_bias(acc, rating, seed=seed)
        if dprime is None:
            dprime, *_ = compute_sdt_resp(stimulus, response)
        n_ratings = int(np.unique(rating).size)
        expected  = _gaussian_meta_I_simple(abs(dprime), n_ratings=n_ratings)
        return float("nan") if expected == 0 else float(mi / expected)

    if backend == "statconfr":
        table = _build_contingency_table(stimulus, response, rating)
        if dprime is None:
            dprime = _estimate_dprime_statconfr(table)
        if bias_correction:
            fn   = lambda t: _sc_meta_Ir1(t, dprime)
            bias = _statconfr_bias(table, fn, seed=seed)
            return float(_sc_meta_Ir1(table, dprime) - bias)
        return _sc_meta_Ir1(table, dprime)

    raise ValueError(f"Unknown backend {backend!r}.")

meta_Ir1_acc(stimulus, response, rating, *, backend='simple', bias_correction=False, seed=42)

Accuracy-normalised relative efficiency (meta-I₁ʳ_acc).

Normalises by Gaussian meta-I at d′ derived from observed accuracy, avoiding SDT assumptions about criterion placement.

Source code in src/metasignal/itmc/measures.py 
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
def meta_Ir1_acc(
    stimulus: np.ndarray,
    response: np.ndarray,
    rating: np.ndarray,
    *,
    backend: Backend = "simple",
    bias_correction: bool = False,
    seed: int = 42,
) -> float:
    """Accuracy-normalised relative efficiency (meta-I₁ʳ_acc).

    Normalises by Gaussian meta-I at d′ derived from observed accuracy,
    avoiding SDT assumptions about criterion placement.
    """
    stimulus = np.asarray(stimulus)
    response = np.asarray(response)
    rating   = np.asarray(rating)

    if backend == "simple":
        acc = (stimulus == response).astype(int)
        mi  = _simple_meta_I(acc, rating)
        if bias_correction:
            mi -= _simple_bias(acc, rating, seed=seed)
        p_acc = float(acc.mean())
        if p_acc <= 0 or p_acc >= 1:
            return float("nan")
        dp_acc    = 2.0 * norm.ppf(p_acc)
        n_ratings = int(np.unique(rating).size)
        expected  = _gaussian_meta_I_simple(abs(dp_acc), n_ratings=n_ratings)
        return float("nan") if expected == 0 else float(mi / expected)

    if backend == "statconfr":
        table = _build_contingency_table(stimulus, response, rating)
        if bias_correction:
            bias = _statconfr_bias(table, _sc_meta_Ir1_acc, seed=seed)
            return float(_sc_meta_Ir1_acc(table) - bias)
        return _sc_meta_Ir1_acc(table)

    raise ValueError(f"Unknown backend {backend!r}.")

meta_Ir2(stimulus, response, rating, *, backend='simple', bias_correction=False, seed=42)

Entropy-normalised absolute metacognitive efficiency (meta-I₂ʳ).

meta-I₂ʳ = meta-I / H₂(accuracy). Range [0, 1].

Both backends use H₂(accuracy) as the denominator; the numerator (meta-I) differs between backends.

Source code in src/metasignal/itmc/measures.py 
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
def meta_Ir2(
    stimulus: np.ndarray,
    response: np.ndarray,
    rating: np.ndarray,
    *,
    backend: Backend = "simple",
    bias_correction: bool = False,
    seed: int = 42,
) -> float:
    """Entropy-normalised absolute metacognitive efficiency (meta-I₂ʳ).

    meta-I₂ʳ = meta-I / H₂(accuracy).  Range [0, 1].

    Both backends use H₂(accuracy) as the denominator; the numerator
    (meta-I) differs between backends.
    """
    stimulus = np.asarray(stimulus)
    response = np.asarray(response)
    rating   = np.asarray(rating)

    if backend == "simple":
        acc = (stimulus == response).astype(int)
        H_r = float(_h2(acc.mean()))
        if H_r == 0:
            return float("nan")
        mi = _simple_meta_I(acc, rating)
        if bias_correction:
            mi -= _simple_bias(acc, rating, seed=seed)
        return float(mi / H_r)

    if backend == "statconfr":
        table = _build_contingency_table(stimulus, response, rating)
        if bias_correction:
            bias = _statconfr_bias(table, _sc_meta_Ir2, seed=seed)
            return float(_sc_meta_Ir2(table) - bias)
        return _sc_meta_Ir2(table)

    raise ValueError(f"Unknown backend {backend!r}.")

RMI(stimulus, response, rating, *, backend='simple', bias_correction=False, seed=42)

Relative Meta-Information, range [0, 1].

RMI = (meta-I − I_min) / (I_max − I_min)

Parameters

backend: 'simple': bounds are H_s − H_r (lower) and H_s (upper) where H_s = H₂(P(S2)) and H_r = H₂(accuracy). These bounds assume I(S;R) as the numerator; since the simple backend uses MI(acc;conf) instead, RMI can fall outside [0, 1] — treat as approximate. Prefer the 'statconfr' backend for RMI.

``'statconfr'``: analytic Fano-inequality bounds
(``get_lower_info_for_one`` / ``get_upper_info_for_one`` in statConfR)
applied to I(stimulus; graded_response).  RMI is properly bounded in
[0, 1] and matches the R package exactly.
Source code in src/metasignal/itmc/measures.py 
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
def RMI(
    stimulus: np.ndarray,
    response: np.ndarray,
    rating: np.ndarray,
    *,
    backend: Backend = "simple",
    bias_correction: bool = False,
    seed: int = 42,
) -> float:
    """Relative Meta-Information, range [0, 1].

    RMI = (meta-I − I_min) / (I_max − I_min)

    Parameters
    ----------
    backend:
        ``'simple'``: bounds are H_s − H_r (lower) and H_s (upper) where
        H_s = H₂(P(S2)) and H_r = H₂(accuracy).  These bounds assume
        I(S;R) as the numerator; since the simple backend uses MI(acc;conf)
        instead, RMI can fall outside [0, 1] — treat as approximate.
        Prefer the ``'statconfr'`` backend for RMI.

        ``'statconfr'``: analytic Fano-inequality bounds
        (``get_lower_info_for_one`` / ``get_upper_info_for_one`` in statConfR)
        applied to I(stimulus; graded_response).  RMI is properly bounded in
        [0, 1] and matches the R package exactly.
    """
    stimulus = np.asarray(stimulus)
    response = np.asarray(response)
    rating   = np.asarray(rating)

    if backend == "simple":
        acc       = (stimulus == response).astype(int)
        p_correct = float(acc.mean())
        p_plus    = float((stimulus == 1).mean())
        H_s       = float(_h2(p_plus))
        H_r       = float(_h2(p_correct))
        min_info  = H_s - H_r
        max_info  = H_s
        denom     = max_info - min_info
        if denom <= 0:
            return float("nan")
        mi = _simple_meta_I(acc, rating)
        if bias_correction:
            mi -= _simple_bias(acc, rating, seed=seed)
        return float((mi - min_info) / denom)

    if backend == "statconfr":
        table = _build_contingency_table(stimulus, response, rating)
        if bias_correction:
            bias = _statconfr_bias(table, _sc_RMI, seed=seed)
            return float(_sc_RMI(table) - bias)
        return _sc_RMI(table)

    raise ValueError(f"Unknown backend {backend!r}.")

permtest_meta_I(stimulus, response, rating, *, backend='simple', n_perm=1000, seed=42)

Test whether meta-I is significantly greater than chance.

Builds a permutation null by shuffling confidence ratings across trials, which breaks the confidence–accuracy link while preserving everything else. Returns the observed value, bias-corrected value (observed − mean(null)), a one-tailed p-value, and null distribution summary statistics.

This is the IT analogue of asking "is metacognition above chance?" without fitting a parametric model.

Parameters

stimulus, response, rating : Trial-level arrays (same conventions as meta_I). backend : 'simple' or 'statconfr'. Both use rating permutation for the null; the backend controls how meta-I is computed on each permuted sample. n_perm : Number of permutations for the null distribution. seed : RNG seed.

Returns

dict with keys: observed — raw meta-I corrected — bias-corrected meta-I (observed − mean(null)) p_value — proportion of null samples ≥ observed (one-tailed) null_mean — mean of null distribution null_std — std of null distribution null — full null array (length n_perm) backend — backend used n_perm — number of permutations used

Source code in src/metasignal/itmc/measures.py 
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
def permtest_meta_I(
    stimulus: np.ndarray,
    response: np.ndarray,
    rating: np.ndarray,
    *,
    backend: Backend = "simple",
    n_perm: int = 1000,
    seed: int = 42,
) -> dict:
    """Test whether meta-I is significantly greater than chance.

    Builds a permutation null by shuffling confidence ratings across trials,
    which breaks the confidence–accuracy link while preserving everything else.
    Returns the observed value, bias-corrected value (observed − mean(null)),
    a one-tailed p-value, and null distribution summary statistics.

    This is the IT analogue of asking "is metacognition above chance?" without
    fitting a parametric model.

    Parameters
    ----------
    stimulus, response, rating :
        Trial-level arrays (same conventions as ``meta_I``).
    backend :
        ``'simple'`` or ``'statconfr'``.  Both use rating permutation for the
        null; the backend controls how meta-I is computed on each permuted sample.
    n_perm :
        Number of permutations for the null distribution.
    seed :
        RNG seed.

    Returns
    -------
    dict with keys:
        ``observed``   — raw meta-I
        ``corrected``  — bias-corrected meta-I (observed − mean(null))
        ``p_value``    — proportion of null samples ≥ observed (one-tailed)
        ``null_mean``  — mean of null distribution
        ``null_std``   — std of null distribution
        ``null``       — full null array (length n_perm)
        ``backend``    — backend used
        ``n_perm``     — number of permutations used
    """
    stimulus = np.asarray(stimulus)
    response = np.asarray(response)
    rating   = np.asarray(rating)
    rng      = np.random.default_rng(seed)

    if backend == "simple":
        acc      = (stimulus == response).astype(int)
        observed = _simple_meta_I(acc, rating)
        null     = np.array([_simple_meta_I(acc, rng.permutation(rating))
                             for _ in range(n_perm)])

    elif backend == "statconfr":
        observed = _sc_meta_I(_build_contingency_table(stimulus, response, rating))
        null_vals = []
        for _ in range(n_perm):
            perm_rating = rng.permutation(rating)
            t = _build_contingency_table(stimulus, response, perm_rating)
            acc_t = _get_accuracy_from_table(t)
            if acc_t is None or np.isnan(acc_t) or abs(acc_t - 1.0) < 1e-6:
                continue
            null_vals.append(_sc_meta_I(t))
        null = np.array(null_vals)

    else:
        raise ValueError(f"Unknown backend {backend!r}. Choose 'simple' or 'statconfr'.")

    null_mean = float(null.mean()) if len(null) else float("nan")
    null_std  = float(null.std())  if len(null) else float("nan")
    p_value   = (int(np.sum(null >= observed)) + 1) / (len(null) + 1) if len(null) else float("nan")
    corrected = observed - null_mean

    return {
        "observed":  float(observed),
        "corrected": float(corrected),
        "p_value":   p_value,
        "null_mean": null_mean,
        "null_std":  null_std,
        "null":      null,
        "backend":   backend,
        "n_perm":    len(null),
    }