Make ltriangle density floor relative to scale

src/TMB/ll_ltriangle.hpp

@@ -85,8 +85,11 @@ Type ll_ltriangle(objective_function<Type>* obj) {
         Type absdev = CppAD::CondExpGe(dev, Type(0.0), dev, -dev);
         Type inside = scalelog - absdev;
         // floor the density argument to keep the likelihood finite (and act as a
-        // soft barrier) if a candidate support fails to cover a data point.
-        Type floored = CppAD::CondExpGt(inside, Type(1e-8), inside, Type(1e-8));
+        // soft barrier) if a candidate support fails to cover a data point. The
+        // floor is scaled to `scalelog` (the peak density argument) so the soft
+        // barrier behaves consistently regardless of the magnitude of the data.
+        Type floor = Type(1e-8) * scalelog;
+        Type floored = CppAD::CondExpGt(inside, floor, inside, floor);
         nll -= weight(i) * (log(floored) - Type(2.0) * log(scalelog) - log(left(i)));
      };
      if(left(i) < right(i)){    // censored values

tests/testthat/test-ltriangle.R

@@ -25,3 +25,34 @@ test_that("ltriangle", {
     expect_snapshot_value(ssd_rltriangle(2), style = "deparse")
   })
 })
+
+test_that("ltriangle fit is invariant to scaling the concentrations", {
+  data <- ssddata::ccme_boron
+  fit <- ssd_fit_dists(data, dists = "ltriangle")
+
+  data_scaled <- data
+  data_scaled$Conc <- data_scaled$Conc * 1000
+  fit_scaled <- ssd_fit_dists(data_scaled, dists = "ltriangle")
+
+  est <- estimates(fit)
+  est_scaled <- estimates(fit_scaled)
+
+  # scaling concentrations by a constant only shifts locationlog by its log and
+  # leaves scalelog unchanged, so hazard concentrations scale by the constant.
+  # The scale-relative density floor keeps this invariance across magnitudes.
+  expect_equal(
+    est_scaled$ltriangle.scalelog,
+    est$ltriangle.scalelog,
+    tolerance = 1e-5
+  )
+  expect_equal(
+    est_scaled$ltriangle.locationlog,
+    est$ltriangle.locationlog + log(1000),
+    tolerance = 1e-5
+  )
+  expect_equal(
+    ssd_hc(fit_scaled, proportion = 0.05)$est,
+    ssd_hc(fit, proportion = 0.05)$est * 1000,
+    tolerance = 1e-4
+  )
+})