diff --git a/crates/mlsirm-core/src/fitstats.rs b/crates/mlsirm-core/src/fitstats.rs
index 9c56e0df9..d6ec0a5ba 100644
--- a/crates/mlsirm-core/src/fitstats.rs
+++ b/crates/mlsirm-core/src/fitstats.rs
@@ -2436,23 +2436,23 @@ pub fn poly_local_dependence(
     model: crate::poly::PolyModel,
     q_theta: usize,
 ) -> Result<PolyLdResult, String> {
-    use crate::poly::{gpcm_logprobs, grm_logprobs, PolyModel};
+    use crate::poly::{gpcm_logprobs, grm_logprobs, validate_poly_item_parameters, PolyModel};
     if n_items < 2 {
         return Err("local dependence needs at least 2 items".into());
     }
     if n_cat < 2 {
         return Err("n_cat must be >= 2".into());
     }
-    if y.len() != n_persons * n_items {
+    let expected_len = crate::checked_mul_usize(
+        n_persons,
+        n_items,
+        "n_persons * n_items exceeds the response buffer size",
+    )?;
+    if y.len() != expected_len {
         return Err("y must have length n_persons * n_items".into());
     }
     validate_optional_observed_length(observed, y.len())?;
-    if slope.len() != n_items {
-        return Err("slope must have length n_items".into());
-    }
-    if cat_params.len() != n_items * (n_cat - 1) {
-        return Err("cat_params must have length n_items*(n_cat-1)".into());
-    }
+    validate_poly_item_parameters(slope, cat_params, n_items, n_cat, model)?;
     validate_observed_categories(y, observed, n_cat)?;
     let z = n_cat - 1;
 
@@ -2603,23 +2603,23 @@ pub fn poly_m2(
     model: crate::poly::PolyModel,
     q_theta: usize,
 ) -> Result<M2Result, String> {
-    use crate::poly::{gpcm_logprobs, grm_logprobs, PolyModel};
+    use crate::poly::{gpcm_logprobs, grm_logprobs, validate_poly_item_parameters, PolyModel};
     if n_items < 3 {
         return Err("M2 needs at least 3 items".into());
     }
     if n_cat < 2 {
         return Err("n_cat must be >= 2".into());
     }
-    if y.len() != n_persons * n_items {
+    let expected_len = crate::checked_mul_usize(
+        n_persons,
+        n_items,
+        "n_persons * n_items exceeds the response buffer size",
+    )?;
+    if y.len() != expected_len {
         return Err("y must have length n_persons * n_items".into());
     }
     validate_optional_observed_length(observed, y.len())?;
-    if slope.len() != n_items {
-        return Err("slope must have length n_items".into());
-    }
-    if cat_params.len() != n_items * (n_cat - 1) {
-        return Err("cat_params must have length n_items*(n_cat-1)".into());
-    }
+    validate_poly_item_parameters(slope, cat_params, n_items, n_cat, model)?;
     validate_observed_categories(y, observed, n_cat)?;
 
     let z = n_cat - 1; // highest threshold index
diff --git a/crates/mlsirm-core/src/poly.rs b/crates/mlsirm-core/src/poly.rs
index b3357be6a..e0c48d9ed 100644
--- a/crates/mlsirm-core/src/poly.rs
+++ b/crates/mlsirm-core/src/poly.rs
@@ -31,6 +31,38 @@ fn validate_observed_categories(
     Ok(())
 }
 
+pub(crate) fn validate_poly_item_parameters(
+    slope: &[f64],
+    cat_params: &[f64],
+    n_items: usize,
+    n_cat: usize,
+    model: PolyModel,
+) -> Result<(), String> {
+    let expected_cat_params = crate::checked_mul_usize(
+        n_items,
+        n_cat - 1,
+        "n_items * (n_cat - 1) exceeds the category-parameter buffer size",
+    )?;
+    if slope.len() != n_items || cat_params.len() != expected_cat_params {
+        return Err("slope/cat_params must match n_items and n_cat".into());
+    }
+    if slope
+        .iter()
+        .chain(cat_params.iter())
+        .any(|value| !value.is_finite())
+    {
+        return Err("slope/cat_params must be finite".into());
+    }
+    if model == PolyModel::Grm {
+        for thresholds in cat_params.chunks_exact(n_cat - 1) {
+            if thresholds.windows(2).any(|pair| pair[0] <= pair[1]) {
+                return Err("GRM thresholds must be strictly decreasing within each item".into());
+            }
+        }
+    }
+    Ok(())
+}
+
 #[inline]
 fn log_sigmoid(x: f64) -> f64 {
     if x >= 0.0 {
@@ -887,12 +919,7 @@ pub fn poly_person_fit(
     if n_cat < 2 {
         return Err("n_cat must be >= 2".into());
     }
-    if slope.len() != n_items {
-        return Err("slope must have length n_items".into());
-    }
-    if cat_params.len() != n_items * (n_cat - 1) {
-        return Err("cat_params must have length n_items*(n_cat-1)".into());
-    }
+    validate_poly_item_parameters(slope, cat_params, n_items, n_cat, model)?;
     if !(prior_sd > 0.0) {
         return Err("prior_sd must be positive".into());
     }
@@ -1036,8 +1063,9 @@ pub fn poly_cat_simulate(
     if n_cat < 2 {
         return Err("n_cat must be >= 2".into());
     }
-    if slope.len() != n_items || cat_params.len() != n_items * (n_cat - 1) {
-        return Err("slope/cat_params must match n_items and n_cat".into());
+    validate_poly_item_parameters(slope, cat_params, n_items, n_cat, model)?;
+    if true_theta.iter().any(|value| !value.is_finite()) {
+        return Err("true_theta must be finite".into());
     }
     if n_items < 2 {
         return Err("CAT needs a bank of at least 2 items".into());
@@ -1885,9 +1913,7 @@ pub fn u3_poly_bootstrap_cutoff(
     if n_cat < 2 {
         return Err("n_cat must be >= 2".into());
     }
-    if slope.len() != n_items || cat_params.len() != n_items * (n_cat - 1) {
-        return Err("slope/cat_params must match n_items and n_cat".into());
-    }
+    validate_poly_item_parameters(slope, cat_params, n_items, n_cat, model)?;
     if n_persons < 1 || n_items < 1 {
         return Err("need at least one person and item".into());
     }
@@ -1918,8 +1944,18 @@ pub fn u3_poly_bootstrap_cutoff(
             PolyModel::Grm => grm_logprobs(base, cp),
         }
     };
-    let mut pool: Vec<f64> = Vec::with_capacity(n_rep * n_persons);
-    let mut y = vec![0usize; n_persons * n_items];
+    let pool_capacity = crate::checked_mul_usize(
+        n_rep,
+        n_persons,
+        "n_rep * n_persons exceeds the bootstrap buffer size",
+    )?;
+    let response_len = crate::checked_mul_usize(
+        n_persons,
+        n_items,
+        "n_persons * n_items exceeds the response buffer size",
+    )?;
+    let mut pool: Vec<f64> = Vec::with_capacity(pool_capacity);
+    let mut y = vec![0usize; response_len];
     for _rep in 0..n_rep {
         for p in 0..n_persons {
             let u1 = u().max(1e-12);
@@ -1942,10 +1978,9 @@ pub fn u3_poly_bootstrap_cutoff(
         let res = u3_poly_person_fit(&y, None, n_persons, n_items, n_cat, None)?;
         pool.extend(res.u3poly.into_iter().filter(|v| v.is_finite()));
     }
-    debug_assert!(
-        !pool.is_empty(),
-        "validated complete bootstrap samples have finite boundary U3 values"
-    );
+    if pool.is_empty() {
+        return Err("bootstrap produced no finite U3 values".into());
+    }
     pool.sort_by(|a, b| a.partial_cmp(b).unwrap());
     let np = pool.len();
     let idx = (np as f64 - 1.0) * (1.0 - alpha);
@@ -2036,16 +2071,9 @@ pub fn poly_information_curves(
     if theta.is_empty() {
         return Err("theta must be non-empty".into());
     }
-    let expected_cat_params =
-        crate::checked_mul_usize(n_items, n_cat - 1, "n_items * (n_cat - 1) overflows usize")?;
-    if slope.len() != n_items || cat_params.len() != expected_cat_params {
-        return Err("slope/cat_params sizes inconsistent with n_items/n_cat".into());
-    }
-    if theta.iter().any(|value| !value.is_finite())
-        || slope.iter().any(|value| !value.is_finite())
-        || cat_params.iter().any(|value| !value.is_finite())
-    {
-        return Err("theta, slope, and cat_params must be finite".into());
+    validate_poly_item_parameters(slope, cat_params, n_items, n_cat, model)?;
+    if theta.iter().any(|value| !value.is_finite()) {
+        return Err("theta must be finite".into());
     }
     let output_len = crate::checked_mul_usize(theta.len(), n_items, "output size overflows")?;
     let mut out = vec![0.0_f64; output_len];
@@ -2094,14 +2122,7 @@ pub fn score_poly_eap(
             return Err("observed must have length n_persons * n_items".into());
         }
     }
-    let n_params =
-        crate::checked_mul_usize(n_items, n_cat - 1, "n_items * (n_cat - 1) overflows usize")?;
-    if slope.len() != n_items || cat_params.len() != n_params {
-        return Err("slope/cat_params sizes inconsistent with n_items/n_cat".into());
-    }
-    if slope.iter().any(|v| !v.is_finite()) || cat_params.iter().any(|v| !v.is_finite()) {
-        return Err("slope and cat_params must be finite".into());
-    }
+    validate_poly_item_parameters(slope, cat_params, n_items, n_cat, model)?;
     for (idx, &yc) in y.iter().enumerate() {
         if observed.map_or(true, |o| o[idx]) && yc >= n_cat {
             return Err(format!(
@@ -2245,17 +2266,7 @@ pub fn poly_s_x2(
     if y.len() != n_cells {
         return Err("y must have length n_persons * n_items".into());
     }
-    if slope.len() != n_items {
-        return Err("slope must have length n_items".into());
-    }
-    let n_item_steps = crate::checked_mul_usize(
-        n_items,
-        n_cat - 1,
-        "n_items * (n_cat - 1) overflows usize",
-    )?;
-    if cat_params.len() != n_item_steps {
-        return Err("cat_params must have length n_items * (n_cat - 1)".into());
-    }
+    validate_poly_item_parameters(slope, cat_params, n_items, n_cat, model)?;
     if let Some(o) = observed {
         if o.len() != n_cells {
             return Err("observed must have length n_persons * n_items".into());
diff --git a/crates/mlsirm-core/src/poly_marginal.rs b/crates/mlsirm-core/src/poly_marginal.rs
index 5d99a674a..fa88e091a 100644
--- a/crates/mlsirm-core/src/poly_marginal.rs
+++ b/crates/mlsirm-core/src/poly_marginal.rs
@@ -2,9 +2,10 @@
 //! EM — the Rust compute path for the GRM/GPCM cell embedded in an interaction
 //! map. Unidimensional trait `theta` with a `latent_dim`-dimensional latent
 //! space; the person latent position `xi` is integrated over a tensor
-//! Gauss-Hermite grid and the item position `zeta_i` is estimated. This is the
-//! `fixed_gamma = 1` identification of Go et al. (2024) lsirm12pl (the distance
-//! weight is fixed to standardize the map scale).
+//! Gauss-Hermite grid and the item position `zeta_i` is estimated. The distance
+//! weight is fixed at 1 as this crate's scale-identification choice; Go et al.'s
+//! lsirm12pl keeps a `gamma` term for continuous responses and leaves ordinal
+//! LSIRM as future work.
 //!
 //! Fully additive: reuses the [`crate::poly`] cells/gradients and the exact
 //! `d eta / d zeta` distance derivative from the binary M-step
@@ -227,11 +228,16 @@ pub fn fit_poly_lsirm(
     if latent_dim < 1 || latent_dim > 3 {
         return Err("latent_dim must be 1..3 for the tensor grid".into());
     }
-    if y.len() != n_persons * n_items {
+    let n_cells = crate::checked_mul_usize(
+        n_persons,
+        n_items,
+        "n_persons * n_items exceeds the response buffer size",
+    )?;
+    if y.len() != n_cells {
         return Err("y must have length n_persons * n_items".into());
     }
     if let Some(o) = observed {
-        if o.len() != n_persons * n_items {
+        if o.len() != n_cells {
             return Err("observed must have length n_persons * n_items".into());
         }
     }
diff --git a/python/fast_mlsirm/polytomous.py b/python/fast_mlsirm/polytomous.py
index 595566287..525b0e4cf 100644
--- a/python/fast_mlsirm/polytomous.py
+++ b/python/fast_mlsirm/polytomous.py
@@ -386,8 +386,9 @@ def fit_lsirm_polytomous(
 ) -> PolyLsirmFit:
     """Fit a latent-space polytomous LSIRM (GRM/GPCM cell in an interaction map)
     by marginal EM — all compute in the Rust core (``poly_marginal``). The
-    distance weight is fixed to 1 (Go et al. 2024 identification); positions are
-    identified up to rotation/reflection/translation. ``NaN`` marks missing.
+    distance weight is fixed to 1 as this crate's scale-identification choice;
+    positions are identified up to rotation/reflection/translation. ``NaN`` marks
+    missing.
     ``n_cat`` is limited to 2..64 and ``max_iter`` to 1..100,000.
     """
     m = str(model).lower()
diff --git a/tests/unit/fitstats_tests.rs b/tests/unit/fitstats_tests.rs
index 91fe15afe..d9139e290 100644
--- a/tests/unit/fitstats_tests.rs
+++ b/tests/unit/fitstats_tests.rs
@@ -876,9 +876,33 @@ fn fitstats_public_boundaries_and_interaction_paths() {
     .is_err());
     assert!(poly_local_dependence(&[0, 0], None, 1, 2, 2, &[1.0], &[0.0, 0.0], pm, 7).is_err());
     assert!(poly_local_dependence(&[0, 0], None, 1, 2, 2, &[1.0, 1.0], &[0.0], pm, 7).is_err());
+    assert!(
+        poly_local_dependence(&[0, 0], None, 1, 2, 2, &[f64::NAN, 1.0], &[0.0, 0.0], pm, 7)
+            .is_err()
+    );
     assert!(
         poly_local_dependence(&[0, 2], None, 1, 2, 2, &[1.0, 1.0], &[0.0, 0.0], pm, 7).is_err()
     );
+    let overflow_ld = std::panic::catch_unwind(|| {
+        poly_local_dependence(
+            &[],
+            None,
+            usize::MAX / 2 + 1,
+            2,
+            2,
+            &[1.0, 1.0],
+            &[0.0, 0.0],
+            pm,
+            7,
+        )
+    });
+    assert!(
+        overflow_ld.is_ok(),
+        "overflowed LD shape must return an error instead of panicking"
+    );
+    // Reads crate Result. Kills the mutation that multiplies n_persons*n_items
+    // directly before checking response length.
+    assert!(overflow_ld.unwrap().is_err());
     let masked_ld = poly_local_dependence(
         &[0, 99],
         Some(&[true, false]),
@@ -915,7 +939,7 @@ fn fitstats_public_boundaries_and_interaction_paths() {
             2,
             3,
             &[1.0, 0.8],
-            &[-0.5, 0.5, -0.25, 0.75],
+            &[0.5, -0.5, 0.75, -0.25],
             model,
             7,
         )
@@ -944,7 +968,28 @@ fn fitstats_public_boundaries_and_interaction_paths() {
     .is_err());
     assert!(poly_m2(&[0, 0, 0], None, 1, 3, 2, &[1.0; 2], &[0.0; 3], pm, 7).is_err());
     assert!(poly_m2(&[0, 0, 0], None, 1, 3, 2, &[1.0; 3], &[0.0; 2], pm, 7).is_err());
+    assert!(poly_m2(&[0, 0, 0], None, 1, 3, 2, &[f64::NAN; 3], &[0.0; 3], pm, 7).is_err());
     assert!(poly_m2(&[0, 0, 2], None, 1, 3, 2, &[1.0; 3], &[0.0; 3], pm, 7).is_err());
+    let overflow_m2 = std::panic::catch_unwind(|| {
+        poly_m2(
+            &[],
+            None,
+            usize::MAX / 3 + 1,
+            3,
+            2,
+            &[1.0; 3],
+            &[0.0; 3],
+            pm,
+            7,
+        )
+    });
+    assert!(
+        overflow_m2.is_ok(),
+        "overflowed M2 shape must return an error instead of panicking"
+    );
+    // Reads crate Result. Kills the mutation that multiplies n_persons*n_items
+    // directly before checking response length.
+    assert!(overflow_m2.unwrap().is_err());
     let mut masked_y: Vec<usize> = (0..20)
         .flat_map(|p| [p % 2, (p / 2) % 2, (p / 3) % 2, (p / 5) % 2])
         .collect();
diff --git a/tests/unit/poly_marginal_tests.rs b/tests/unit/poly_marginal_tests.rs
index af86ce692..55a52cfb3 100644
--- a/tests/unit/poly_marginal_tests.rs
+++ b/tests/unit/poly_marginal_tests.rs
@@ -75,6 +75,28 @@ fn poly_marginal_boundaries_and_grm_paths_are_explicit() {
     assert!(fit_poly_lsirm(&[], None, 0, 1, 2, 1, PolyModel::Grm, 7, 7, 1, 1e-6).is_err());
     assert!(fit_poly_lsirm(&[], None, 1, 0, 2, 1, PolyModel::Grm, 7, 7, 1, 1e-6).is_err());
     assert!(fit_poly_lsirm(&[], None, 1, 1, 2, 1, PolyModel::Grm, 7, 7, 1, 1e-6).is_err());
+    let overflow_lsirm = std::panic::catch_unwind(|| {
+        fit_poly_lsirm(
+            &[],
+            None,
+            usize::MAX / 2 + 1,
+            2,
+            2,
+            1,
+            PolyModel::Grm,
+            7,
+            7,
+            1,
+            1e-6,
+        )
+    });
+    assert!(
+        overflow_lsirm.is_ok(),
+        "overflowed poly LSIRM shape must return an error instead of panicking"
+    );
+    // Reads crate Result. Kills the mutation that multiplies n_persons*n_items
+    // directly before checking response length.
+    assert!(overflow_lsirm.unwrap().is_err());
     assert!(fit_poly_lsirm(&[0], Some(&[]), 1, 1, 2, 1, PolyModel::Grm, 7, 7, 1, 1e-6,).is_err());
     assert!(fit_poly_lsirm(
         &[3],
diff --git a/tests/unit/poly_tests.rs b/tests/unit/poly_tests.rs
index 986f8d36d..1aca59fdd 100644
--- a/tests/unit/poly_tests.rs
+++ b/tests/unit/poly_tests.rs
@@ -259,6 +259,44 @@ fn poly_public_boundaries_and_small_diagnostic_paths() {
     ] {
         assert!(result.is_err());
     }
+    let overflow_person_fit = std::panic::catch_unwind(|| {
+        poly_person_fit(
+            &[],
+            None,
+            1,
+            2,
+            usize::MAX / 2 + 2,
+            &[1.0, 1.0],
+            &[],
+            PolyModel::Gpcm,
+            7,
+            0.0,
+            1.0,
+            -2.0,
+        )
+    });
+    assert!(
+        overflow_person_fit.is_ok(),
+        "overflowed poly person-fit item shape must return an error instead of panicking"
+    );
+    // Reads crate Result. Kills the mutation that multiplies n_items*(n_cat-1)
+    // directly before validating category-parameter length.
+    assert!(overflow_person_fit.unwrap().is_err());
+    let equal_grm_thresholds = score_poly_eap(
+        &[1],
+        None,
+        1,
+        1,
+        3,
+        &[1.0],
+        &[0.0, 0.0],
+        PolyModel::Grm,
+        7,
+    );
+    let equal_grm_error = equal_grm_thresholds.expect_err("equal GRM thresholds accepted");
+    // Reads crate error text. Kills the mutation that permits a zero-probability
+    // middle category and lets Result APIs return Ok(NaN).
+    assert!(equal_grm_error.contains("thresholds"));
     let sparse_observed = [false, false, true, true];
     let person_fit = poly_person_fit(
         &y,
@@ -597,6 +635,66 @@ fn poly_public_boundaries_and_small_diagnostic_paths() {
     ] {
         assert!(result.is_err());
     }
+    let overflow_u3 =
+        std::panic::catch_unwind(|| u3_poly_person_fit(&[], None, usize::MAX / 2 + 1, 2, 3, None));
+    assert!(
+        overflow_u3.is_ok(),
+        "overflowed U3 shape must return an error instead of panicking"
+    );
+    assert!(overflow_u3.unwrap().is_err());
+    let empty_bootstrap = std::panic::catch_unwind(|| {
+        u3_poly_bootstrap_cutoff(1, 1, 3, &[0.0], &[0.0, 0.0], PolyModel::Gpcm, 0.05, 1, 1)
+    });
+    assert!(
+        empty_bootstrap.is_ok(),
+        "U3 bootstrap must return an error when all simulated U3 values are non-finite"
+    );
+    assert!(empty_bootstrap.unwrap().is_err());
+    let nan_cutoff = u3_poly_bootstrap_cutoff(
+        1,
+        1,
+        3,
+        &[f64::NAN],
+        &[0.0, 0.0],
+        PolyModel::Gpcm,
+        0.05,
+        1,
+        7,
+    );
+    assert!(nan_cutoff.is_err());
+    let unordered_grm_cutoff = u3_poly_bootstrap_cutoff(
+        1,
+        1,
+        3,
+        &[1.0],
+        &[-1.0, 1.0],
+        PolyModel::Grm,
+        0.05,
+        1,
+        7,
+    );
+    let unordered_grm_error = unordered_grm_cutoff.expect_err("unordered GRM thresholds accepted");
+    // Reads crate error text. Kills the mutation that bypasses shared GRM
+    // threshold validation and lets grm_logprobs produce non-finite bootstrap draws.
+    assert!(unordered_grm_error.contains("thresholds"));
+    let overflow_cutoff = std::panic::catch_unwind(|| {
+        u3_poly_bootstrap_cutoff(
+            usize::MAX / 2 + 1,
+            2,
+            3,
+            &[0.0, 0.0],
+            &[0.0; 4],
+            PolyModel::Gpcm,
+            0.05,
+            1,
+            7,
+        )
+    });
+    assert!(
+        overflow_cutoff.is_ok(),
+        "U3 cutoff must return an error for overflowed bootstrap buffers"
+    );
+    assert!(overflow_cutoff.unwrap().is_err());
     let none_observed = [false, false, true, false];
     let y_masked = [99usize, 99, 1, 99];
     let u3 = u3_poly_person_fit(&y_masked, Some(&none_observed), 2, 2, 3, Some(-1.0)).unwrap();
@@ -627,6 +725,8 @@ fn poly_public_boundaries_and_small_diagnostic_paths() {
 
     assert!(poly_information_curves(&[0.0], &slope, &cat, 2, 1, PolyModel::Gpcm).is_err());
     assert!(poly_information_curves(&[0.0], &slope, &[0.0], 2, 3, PolyModel::Gpcm).is_err());
+    let unordered_grm = [-1.0, 1.0, 0.5, -0.5];
+    assert!(poly_information_curves(&[0.0], &slope, &unordered_grm, 2, 3, PolyModel::Grm).is_err());
     for model in [PolyModel::Gpcm, PolyModel::Grm] {
         let information =
             poly_information_curves(&[-1.0, 0.0, 1.0], &slope, &cat, 2, 3, model).unwrap();
@@ -651,6 +751,7 @@ fn poly_public_boundaries_and_small_diagnostic_paths() {
             7,
         ),
         score_poly_eap(&y, None, 2, 2, 3, &[1.0], &cat, PolyModel::Gpcm, 7),
+        score_poly_eap(&y, None, 2, 2, 3, &slope, &unordered_grm, PolyModel::Grm, 7),
     ] {
         assert!(result.is_err());
     }
@@ -736,6 +837,38 @@ fn poly_public_boundaries_and_small_diagnostic_paths() {
     ] {
         assert!(result.is_err());
     }
+    let overflow_sx2 = std::panic::catch_unwind(|| {
+        poly_s_x2(
+            &[],
+            None,
+            usize::MAX / 2 + 1,
+            2,
+            3,
+            &slope,
+            &cat,
+            PolyModel::Gpcm,
+            7,
+            1.0,
+        )
+    });
+    assert!(
+        overflow_sx2.is_ok(),
+        "overflowed poly S-X2 shape must return an error instead of panicking"
+    );
+    assert!(overflow_sx2.unwrap().is_err());
+    assert!(poly_s_x2(
+        &y,
+        None,
+        2,
+        2,
+        3,
+        &[f64::NAN, 1.0],
+        &cat,
+        PolyModel::Gpcm,
+        7,
+        1.0,
+    )
+    .is_err());
     let empty_sx2 = poly_s_x2(
         &[99, 99, 1, 99],
         Some(&[false; 4]),
@@ -2167,6 +2300,63 @@ fn poly_cat_recovers_and_beats_random() {
             (-2.0 * u1.ln()).sqrt() * (2.0 * std::f64::consts::PI * u2).cos()
         })
         .collect();
+    let bad_cat = std::panic::catch_unwind(|| {
+        poly_cat_simulate(
+            &[0.0],
+            &[f64::NAN; 2],
+            &[0.0; 4],
+            2,
+            k,
+            PolyModel::Gpcm,
+            21,
+            0.30,
+            1,
+            2,
+            true,
+            111,
+        )
+    });
+    assert!(
+        bad_cat.is_ok(),
+        "poly_cat_simulate must return an error for non-finite item parameters"
+    );
+    assert!(bad_cat.unwrap().is_err());
+    let bad_grm_cat = std::panic::catch_unwind(|| {
+        poly_cat_simulate(
+            &[0.0],
+            &[1.0; 2],
+            &[-1.0, 1.0, 0.5, -0.5],
+            2,
+            k,
+            PolyModel::Grm,
+            21,
+            0.30,
+            1,
+            2,
+            true,
+            113,
+        )
+    });
+    assert!(
+        bad_grm_cat.is_ok(),
+        "poly_cat_simulate must return an error for unordered GRM thresholds"
+    );
+    assert!(bad_grm_cat.unwrap().is_err());
+    assert!(poly_cat_simulate(
+        &[f64::NAN],
+        &slope,
+        &cat,
+        n_items,
+        k,
+        PolyModel::Gpcm,
+        21,
+        0.30,
+        5,
+        30,
+        false,
+        112,
+    )
+    .is_err());
     // adaptive, variable length: stop at SE < 0.30
     let var = poly_cat_simulate(
         &true_theta,