Point estimation: trustworthy find_MAP, new find_MLE

docs/api/hssm.md

@@ -17,6 +17,9 @@ Use `hssm.HSSM` class to construct an HSSM model.
             - sample_prior_predictive
             - vi
             - find_MAP
+            - find_MLE
+            - map
+            - mle
             - log_likelihood
             - graph
             - plot_predictive

docs/api/point_estimate.md

@@ -0,0 +1,14 @@
+`hssm.PointEstimate` is what [`HSSM.find_MAP`](hssm.md) and
+[`HSSM.find_MLE`](hssm.md) return, and what the `model.map` and `model.mle`
+properties hold. It is a `dict` subclass, so it can be passed anywhere a plain
+point dictionary was accepted before — notably `model.sample(initvals=...)` —
+while also carrying the optimizer metadata and the ArviZ-friendly exporters
+documented below. See the
+[point estimation tutorial](../tutorials/map_mle.ipynb) for worked examples.
+
+::: hssm.PointEstimate
+    options:
+      # The project default suppresses attribute docs, but for this class the
+      # attributes *are* the API — everything the optimizer reports back lives
+      # there rather than on a method.
+      show_docstring_attributes: true

docs/api/rl.md

@@ -25,6 +25,9 @@ sampling model from a named model string.
             - sample_prior_predictive
             - vi
             - find_MAP
+            - find_MLE
+            - map
+            - mle
             - log_likelihood
             - plot_predictive
             - plot_quantile_probability

docs/changelog.md

@@ -6,7 +6,9 @@
 
 2. **`plot_model_cartoon` correctness follow-up** (#1125). Each posterior draw's geometry now derives from **one coherent θ vector** — the across-trials mean by default, or one trial via the new `obs=` parameter (which then conditions every simulated layer on that trial; the observed histogram stays pooled, with a warning). New `random_state=` (int or `np.random.Generator`) makes the whole figure reproducible: which posterior draws are displayed (previously chosen by the unseeded global RNG before any seed applied), every simulator seed, and the trajectories — facets consume successive segments of one stream. Geometry and trajectory simulations now receive a `max_t` horizon derived from the x-limits instead of the simulator's 20 s default (previously ~75% of every boundary polyline lay outside the axes); the noisy RT simulations deliberately keep the long horizon so defective densities are not silently re-normalized by censoring. Trajectories are noisy realizations of the reduced reference θ, so their crossing markers land on the drawn boundary by construction. `n_trajectories`, `xlims`, and `ylims` are promoted to documented parameters (None keeps each renderer's defaults). Bug fix: the >2-choice renderer simulated each draw's RT histogram with `n_samples=1` (silently ignoring `n_reps`) and one seed shared across draws, so its uncertainty bands summarized degenerate, noise-correlated histograms. Deliberate value changes for a fixed `random_state`: a new seeding protocol, trial-mean instead of trial-0 geometry for regression models (intercept-only models unaffected), and reference-θ trajectories. Documented caveat: geometry is nonlinear in θ, so the trial-mean curve is not the mean of per-trial curves and need not sit mid-band — the reduction is a display convention.
 
-3. **`plot_model_cartoon` accepts `hist_height="auto"`** (#1127): fits the tallest RT-histogram curve to 90% of the vertical headroom between the histogram baseline (ribbon-aware — expanding bounds are cleared too) and the upper y-limit, so histograms can never overrun the axes; resolves per facet. Its mirror image, `ylims="auto"`, keeps the raw density scale and grows the frame around the content instead (never below the default limits; the two spellings are mutually exclusive). The defaults remain unchanged. The >2-choice renderer now also draws per-choice drift-uncertainty cones in band mode (graded fills in each accumulator's color; previously slope uncertainty appeared only as per-draw spaghetti in the samples display).
+3. **Point estimation: `find_MAP` is now trustworthy, and `find_MLE` is new** (#1102). `find_MAP` previously delegated to `pm.find_MAP` with no starting point, so PyMC rebuilt its own (`t=2.0`, `a=2.0`) and ignored HSSM's processed initial values. For many models — including a plain hierarchical DDM on `cavanagh_theta` in the default float64 configuration — the *gradient* of the log-density is non-finite there while the logp is finite, so PyMC's start check passed, L-BFGS-B aborted at `nit=0`, and the untouched start point was stored as `model.map` with no warning. `find_MAP` now defaults `start` to `model.initvals`, audits the optimizer's result and **never caches a failed estimate** (it warns and returns `None`, or raises under `strict=True`); `sample(initvals="map")` raises instead of silently falling back when the MAP does not converge. It also refreshes RL/aDDM extra fields the way `sample()` does, warns and switches to the gradient-free `"Powell"` method under `hssm.set_floatX("float32")` (where L-BFGS-B stops early *and reports success*), and gains `n_starts`, `strict`, `se`, `seed` and `method` parameters. New `find_MLE()` maximizes the observed log-likelihood only; it raises on hierarchical models (the group-level scale is unidentified without the priors under the non-centered parameterization — use `find_MAP`) and on models with `pm.Potential` terms, with an `allow_unidentified=True` escape hatch. Both return a `PointEstimate` — a `dict` subclass, so existing `initvals=` code keeps working — carrying `.params`, `.success`, `.logp`, `.se`, `.to_dataframe()` and `.to_datatree()` (a 1-chain/1-draw posterior that feeds `az.summary`, `sample_posterior_predictive` and the plotting functions; `sd` is 0 and `ess`/`r_hat` are `NaN` by construction). New `model.mle` property, `hssm.PointEstimate` exported at the top level (with a `.copy()` that keeps the metadata a plain `dict.copy()` would drop), and `model.map`'s guard changed from a falsy check to `is None`. Both methods take `progressbar=`, and both record the optimizer that *actually* ran — a likelihood without a gradient is silently switched to `"Powell"` by SciPy's caller, so the requested method would otherwise be misreported. `sample(sampler="laplace")` now warns: bambi's laplace path provides no initval plumbing, so it hits the same bad start and HSSM cannot reach it. See the new "MAP and MLE point estimation" tutorial.
+
+4. **`plot_model_cartoon` accepts `hist_height="auto"`** (#1127): fits the tallest RT-histogram curve to 90% of the vertical headroom between the histogram baseline (ribbon-aware — expanding bounds are cleared too) and the upper y-limit, so histograms can never overrun the axes; resolves per facet. Its mirror image, `ylims="auto"`, keeps the raw density scale and grows the frame around the content instead (never below the default limits; the two spellings are mutually exclusive). The defaults remain unchanged. The >2-choice renderer now also draws per-choice drift-uncertainty cones in band mode (graded fills in each accumulator's color; previously slope uncertainty appeared only as per-draw spaghetti in the samples display).
 
 ### 0.4.0