diff --git a/CHANGELOG.md b/CHANGELOG.md
index ab2c9ffb..870098a2 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -5,6 +5,21 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## [0.26.0.dev3] - 2026-10-14
+
+### Added
+
+- Added `benchmarks/hierarchical_profile.py`, which shows where the time of a `VariationalBayes` fit goes and counts the objective evaluations per participant and iteration; against cbm (RLRW, 100 participants, 120 trials), cpm makes about as many evaluations but each one costs more, at 17.3 against 5.0 ms per participant and iteration
+- Added a `hessian` argument to `cpm.optimisation.FminBound`, `Fmin`, `Minimize` and `Bads`; `hessian="finite_differences"` computes the Hessian at the optimum with bounded finite differences (`cpm.core.optimisers.finite_difference_hessian`), cutting Hessian evaluations from 22 to 7 per participant and iteration in the speed comparison. The default, `"numdifftools"`, keeps the previous Hessian
+- Added a `start` argument to `cpm.hierarchical.EmpiricalBayes` and `cpm.hierarchical.VariationalBayes`; `start="prior_mean"` starts the first optimisation of every participant from the current prior's location, clipped into the bounds, as cbm does, so fits stop sooner. The default, `"random"`, keeps the previous starts
+- Added a `positive_definite` column to the `fit` of `cpm.hierarchical.EmpiricalBayes` and `cpm.hierarchical.VariationalBayes`, recording whether each participant's Hessian was positive definite, and a `not_positive_definite` column to their `hyperparameters` with the count; an estimate on a bound can have a non-positive-definite Hessian and is kept, as the bounds require
+- Added the results of the speed comparison with cbm (100 participants, 120 trials) in `benchmarks/results/hierarchical-comparison-*.csv`: combining `hessian="finite_differences"`, `start="prior_mean"` and `number_of_starts=1` brings `EmpiricalBayes`/`VariationalBayes` from 20.8/20.1 to 6.8/7.0 ms per participant and iteration, against cbm's 4.4 ms, with RMSE against 0.26.0.dev1 staying within bootstrap 95% intervals throughout
+
+### Changed
+
+- `cpm.optimisation.FminBound` with `approx_grad=True` now computes the finite-difference gradient itself, together with the objective, exactly as SciPy does, with identical results; with `LogLikelihood.bernoulli`, this makes the hierarchical fits of the speed comparison 19% faster
+- The defaults of the new arguments, `hessian="numdifftools"` and `start="random"`, give the same results as before for the same seed; the faster options have to be chosen
+
 ## [0.26.0.dev2] - 2026-10-06
 
 ### Added
diff --git a/benchmarks/hierarchical_profile.py b/benchmarks/hierarchical_profile.py
new file mode 100644
index 00000000..cc611445
--- /dev/null
+++ b/benchmarks/hierarchical_profile.py
@@ -0,0 +1,261 @@
+"""
+Where does the time of a hierarchical fit go?
+
+    python benchmarks/hierarchical_profile.py [--data DATASET.csv] [--iterations 5] [--label NAME]
+
+Fits one dataset of the four-armed bandit (100 participants, 120 trials) with
+`RLRW` and `VariationalBayes` for a fixed number of iterations, one chain and
+two optimiser starts per participant, under `cProfile`. Prints the share of
+time in the optimiser, the Hessian, `evaluate_fit`, the population update and
+the parts of `objective`, and the number of objective evaluations per
+participant and iteration, separately for the optimiser, the Hessian and
+`evaluate_fit`. Writes the table to benchmarks/results/<label>.csv.
+
+`--data` is a dataset in the format `RLRW` reads (columns ppt, trial, arm_left,
+arm_right, reward_left, reward_right, response). Without it, the task and the
+responses are simulated from a known population. The settings match the speed
+comparison with cbm: bounds alpha [0.01, 1] and temperature [0, 10], a starting
+prior of mean 0.5 and sd 0.5 for alpha and mean 5 and sd 5 for temperature.
+
+The profiler adds overhead to Python-heavy code, so the shares are indicative;
+the fit is also timed once without it.
+"""
+
+import argparse
+import cProfile
+import pathlib
+import pstats
+import sys
+import time
+import warnings
+
+import numpy as np
+import pandas as pd
+from scipy.stats import truncnorm
+
+sys.path.insert(0, str(pathlib.Path(__file__).parent))
+
+import cpm  # noqa: E402
+import cpm.optimisation.fmin as fmin_module  # noqa: E402
+from cpm.applications.reinforcement_learning import RLRW  # noqa: E402
+from cpm.hierarchical import VariationalBayes  # noqa: E402
+from cpm.optimisation import FminBound, minimise  # noqa: E402
+
+RESULTS = pathlib.Path(__file__).parent / "results"
+PARTICIPANTS, TRIALS, STIMULI = 100, 120, 4
+BOUNDS = {"alpha": (1e-2, 1.0), "temperature": (0.0, 10.0)}
+STARTING_PRIOR = {"alpha": {"mean": 0.5, "sd": 0.5}, "temperature": {"mean": 5.0, "sd": 5.0}}
+POPULATION = {"alpha": {"mean": 0.5, "sd": 0.15}, "temperature": {"mean": 4.5, "sd": 1.5}}
+
+
+def simulate(seed=2026):
+    """A random four-armed bandit task with responses of `RLRW` from a known population."""
+    rng = np.random.default_rng(seed)
+    n, trials = PARTICIPANTS, TRIALS
+    arms = np.array([rng.choice(STIMULI, size=2, replace=False) for _ in range(n * trials)]).reshape(n, trials, 2)
+    reward_probability = np.array([0.2, 0.4, 0.6, 0.8])
+    rewards = (rng.random((n, trials, 2)) < reward_probability[arms]).astype(float)
+    parameters = {
+        name: truncnorm.rvs(
+            (BOUNDS[name][0] - p["mean"]) / p["sd"], (BOUNDS[name][1] - p["mean"]) / p["sd"],
+            loc=p["mean"], scale=p["sd"], size=n, random_state=rng,
+        )
+        for name, p in POPULATION.items()
+    }
+    values = np.full((n, STIMULI), 1 / STIMULI)
+    rows = np.arange(n)
+    response = np.empty((n, trials), dtype=int)
+    for t in range(trials):
+        scaled = values[rows[:, None], arms[:, t]] * parameters["temperature"][:, None]
+        scaled -= scaled.max(axis=1, keepdims=True)
+        policy = np.exp(scaled)
+        policy /= policy.sum(axis=1, keepdims=True)
+        choice = (rng.random(n) >= policy[:, 0]).astype(int)
+        stimulus = arms[rows, t, choice]
+        values[rows, stimulus] += parameters["alpha"] * (rewards[rows, t, choice] - values[rows, stimulus])
+        response[:, t] = choice
+    return pd.DataFrame(
+        {
+            "ppt": np.repeat(np.arange(1, n + 1), trials),
+            "trial": np.tile(np.arange(1, trials + 1), n),
+            "arm_left": arms[:, :, 0].ravel() + 1,
+            "arm_right": arms[:, :, 1].ravel() + 1,
+            "reward_left": rewards[:, :, 0].ravel(),
+            "reward_right": rewards[:, :, 1].ravel(),
+            "response": response.ravel(),
+            "observed": response.ravel(),
+        }
+    )
+
+
+def load(path):
+    data = pd.read_csv(path)
+    ppts = np.sort(data.ppt.unique())[:PARTICIPANTS]
+    data = data[data.ppt.isin(ppts) & (data.trial <= TRIALS)].reset_index(drop=True)
+    data["observed"] = data["response"]
+    return data
+
+
+def setup(data, iterations, start="random", **options):
+    """`VariationalBayes` with the settings of the speed comparison, forced to run `iterations` iterations; `options` go to `FminBound`."""
+    first = data[data.ppt == data.ppt.iloc[0]]
+    model = RLRW(data=first, dimensions=STIMULI, parameters_settings=[[0.5, *BOUNDS["alpha"]], [2.0, *BOUNDS["temperature"]]])
+    model.parameters.update_prior(**STARTING_PRIOR)
+    optimiser = FminBound(
+        model=model,
+        data=data.groupby("ppt"),
+        minimisation=minimise.LogLikelihood.bernoulli,
+        prior=True,
+        number_of_starts=2,
+        ppt_identifier="ppt",
+        parallel=False,
+        display=False,
+        approx_grad=True,
+        **options,
+    )
+    return VariationalBayes(
+        optimiser=optimiser,
+        iteration=iterations,
+        chain=1,
+        tolerance_lme=0.0,
+        tolerance_param=0.0,
+        hyperpriors={"a0": np.array([p["mean"] for p in STARTING_PRIOR.values()]), "b": 1, "v": 0.5, "s": np.repeat(0.01, 2)},
+        quiet=True,
+        start=start,
+    )
+
+
+class Counter:
+    """Counts the objective evaluations of `FminBound` by what asked for them."""
+
+    def __init__(self):
+        self.phase = "optimiser"
+        self.counts = {"optimiser": 0, "hessian": 0, "evaluate_fit": 0}
+        self.originals = (fmin_module.objective, fmin_module.numerical_hessian, fmin_module.finite_difference_hessian,
+                          fmin_module.evaluate_fit)
+
+    def __enter__(self):
+        objective, numerical_hessian, finite_difference_hessian, evaluate_fit = self.originals
+
+        def counted_objective(*args, **kwargs):
+            self.counts[self.phase] += 1
+            return objective(*args, **kwargs)
+
+        def counted(hessian):
+            def counted_hessian(*args, **kwargs):
+                self.phase = "hessian"
+                try:
+                    return hessian(*args, **kwargs)
+                finally:
+                    self.phase = "optimiser"
+
+            return counted_hessian
+
+        def counted_evaluate_fit(*args, **kwargs):
+            self.counts["evaluate_fit"] += 1  # one model run per call
+            return evaluate_fit(*args, **kwargs)
+
+        fmin_module.objective = counted_objective
+        fmin_module.numerical_hessian = counted(numerical_hessian)
+        fmin_module.finite_difference_hessian = counted(finite_difference_hessian)
+        fmin_module.evaluate_fit = counted_evaluate_fit
+        return self
+
+    def __exit__(self, *exc):
+        (fmin_module.objective, fmin_module.numerical_hessian, fmin_module.finite_difference_hessian,
+         fmin_module.evaluate_fit) = self.originals
+
+
+def profile(vb):
+    profiler = cProfile.Profile()
+    profiler.enable()
+    vb.optimise()
+    profiler.disable()
+    stats = pstats.Stats(profiler)
+    stats.calc_callees()
+    return stats
+
+
+def find(stats, filename, function):
+    """The (file, line, name) key of a profiled function, by the basename of its file."""
+    for key in stats.stats:
+        if key[2] == function and pathlib.Path(key[0]).name == filename:
+            return key
+    return None
+
+
+def inclusive(stats, filename, function):
+    key = find(stats, filename, function)
+    return stats.stats[key][3] if key else 0.0
+
+
+def inclusive_from(stats, caller, filename, function):
+    """Seconds spent in `function` when called from `caller`."""
+    callees = stats.all_callees.get(caller, {})
+    return sum(value[3] for key, value in callees.items() if key[2] == function and pathlib.Path(key[0]).name == filename)
+
+
+def main():
+    parser = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
+    parser.add_argument("--data", help="a dataset CSV; simulated if not given")
+    parser.add_argument("--iterations", type=int, default=5)
+    parser.add_argument("--label", default="hierarchical-profile")
+    parser.add_argument("--hessian", default="numdifftools", help="the `hessian` option of FminBound")
+    parser.add_argument("--start", default="random", help="the `start` option of VariationalBayes")
+    args = parser.parse_args()
+
+    warnings.simplefilter("ignore")
+    np.seterr(all="ignore")
+    data = load(args.data) if args.data else simulate()
+    participants, iterations = data.ppt.nunique(), args.iterations
+    options = {"hessian": args.hessian, "start": args.start}
+    print(f"cpm {cpm.__version__}: {participants} participants, {data.trial.max()} trials, {iterations} iterations, 2 starts, {options}")
+
+    # the compiled model and everything else that loads on the first call
+    setup(data[data.ppt.isin(data.ppt.unique()[:5])], 1, **options).optimise()
+
+    # once without the profiler, for the time
+    np.random.seed(2026)
+    start = time.perf_counter()
+    with Counter() as counter:
+        setup(data, iterations, **options).optimise()
+    seconds = time.perf_counter() - start
+    per = participants * iterations
+    counts = {k: v / per for k, v in counter.counts.items()}
+
+    np.random.seed(2026)
+    stats = profile(setup(data, iterations, **options))
+    total = inclusive(stats, "variational.py", "optimise")
+    objective_key = find(stats, "optimisers.py", "objective")
+    rows = [
+        ("fit", "seconds without profiler", seconds),
+        ("fit", "ms per participant and iteration", seconds * 1e3 / per),
+        ("evaluations per participant and iteration", "optimiser", counts["optimiser"]),
+        ("evaluations per participant and iteration", "hessian", counts["hessian"]),
+        ("evaluations per participant and iteration", "evaluate_fit", counts["evaluate_fit"]),
+        ("share of profiled time", "fmin_l_bfgs_b", inclusive(stats, "_lbfgsb_py.py", "fmin_l_bfgs_b") / total),
+        ("share of profiled time", "hessian", (inclusive(stats, "optimisers.py", "numerical_hessian")
+                                               + inclusive(stats, "optimisers.py", "finite_difference_hessian")) / total),
+        ("share of profiled time", "evaluate_fit", inclusive(stats, "optimisers.py", "evaluate_fit") / total),
+        ("share of profiled time", "update_population", inclusive(stats, "variational.py", "update_population") / total),
+        ("share of profiled time", "get_lme", inclusive(stats, "variational.py", "get_lme") / total),
+        ("share of profiled time", "objective (all callers)", inclusive(stats, "optimisers.py", "objective") / total),
+        ("share of objective", "SessionWrapper.run", inclusive_from(stats, objective_key, "session.py", "run")),
+        ("share of objective", "  of which the RLRW model function", inclusive(stats, "reinforcement_learning.py", "__call__")),
+        ("share of objective", "Wrapper.reset", inclusive_from(stats, objective_key, "session.py", "reset")),
+        ("share of objective", "deepcopy", inclusive_from(stats, objective_key, "copy.py", "deepcopy")),
+        ("share of objective", "loss", inclusive_from(stats, objective_key, "minimise.py", "bernoulli")),
+        ("share of objective", "Parameters.PDF", inclusive_from(stats, objective_key, "parameters.py", "PDF")),
+    ]
+    objective_seconds = inclusive(stats, "optimisers.py", "objective")
+    rows = [(s, i, v / objective_seconds if s == "share of objective" else v) for s, i, v in rows]
+    table = pd.DataFrame(rows, columns=["section", "item", "value"])
+    table["value"] = table["value"].round(4)
+    print(table.to_string(index=False))
+
+    RESULTS.mkdir(exist_ok=True)
+    table.to_csv(RESULTS / f"{args.label}.csv", index=False)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/benchmarks/results/hierarchical-comparison-individual.csv b/benchmarks/results/hierarchical-comparison-individual.csv
new file mode 100644
index 00000000..99ff9a51
--- /dev/null
+++ b/benchmarks/results/hierarchical-comparison-individual.csv
@@ -0,0 +1,31 @@
+config,method,parameter,individual_rmse
+0.26.0.dev1 (Phase 0 commit),cpm-empirical,alpha,"0.085 [0.073, 0.097]"
+0.26.0.dev1 (Phase 0 commit),cpm-empirical,temperature,"0.763 [0.640, 0.895]"
+0.26.0.dev1 (Phase 0 commit),cpm-variational,alpha,"0.085 [0.074, 0.097]"
+0.26.0.dev1 (Phase 0 commit),cpm-variational,temperature,"0.763 [0.642, 0.895]"
+"new, defaults",cpm-empirical,alpha,"0.085 [0.073, 0.097]"
+"new, defaults",cpm-empirical,temperature,"0.763 [0.640, 0.895]"
+"new, defaults",cpm-variational,alpha,"0.085 [0.074, 0.097]"
+"new, defaults",cpm-variational,temperature,"0.763 [0.642, 0.895]"
+"new, hessian",cpm-empirical,alpha,"0.084 [0.073, 0.096]"
+"new, hessian",cpm-empirical,temperature,"0.760 [0.637, 0.892]"
+"new, hessian",cpm-variational,alpha,"0.085 [0.074, 0.096]"
+"new, hessian",cpm-variational,temperature,"0.758 [0.636, 0.889]"
+"new, start",cpm-empirical,alpha,"0.085 [0.074, 0.097]"
+"new, start",cpm-empirical,temperature,"0.760 [0.637, 0.893]"
+"new, start",cpm-variational,alpha,"0.085 [0.074, 0.097]"
+"new, start",cpm-variational,temperature,"0.762 [0.640, 0.892]"
+"new, start + hessian",cpm-empirical,alpha,"0.084 [0.073, 0.096]"
+"new, start + hessian",cpm-empirical,temperature,"0.759 [0.635, 0.889]"
+"new, start + hessian",cpm-variational,alpha,"0.085 [0.074, 0.096]"
+"new, start + hessian",cpm-variational,temperature,"0.755 [0.635, 0.884]"
+"new, start + hessian, 1 start",cpm-empirical,alpha,"0.085 [0.073, 0.096]"
+"new, start + hessian, 1 start",cpm-empirical,temperature,"0.758 [0.635, 0.890]"
+"new, start + hessian, 1 start",cpm-variational,alpha,"0.085 [0.074, 0.097]"
+"new, start + hessian, 1 start",cpm-variational,temperature,"0.756 [0.634, 0.886]"
+comparison (cpm-empirical),cpm-empirical,alpha,"0.085 [0.073, 0.097]"
+comparison (cpm-empirical),cpm-empirical,temperature,"0.763 [0.639, 0.895]"
+comparison (cpm-variational),cpm-variational,alpha,"0.085 [0.074, 0.097]"
+comparison (cpm-variational),cpm-variational,temperature,"0.764 [0.642, 0.895]"
+comparison (cbm),cbm,alpha,"0.085 [0.074, 0.097]"
+comparison (cbm),cbm,temperature,"0.765 [0.634, 0.906]"
diff --git a/benchmarks/results/hierarchical-comparison-population.csv b/benchmarks/results/hierarchical-comparison-population.csv
new file mode 100644
index 00000000..c61beb33
--- /dev/null
+++ b/benchmarks/results/hierarchical-comparison-population.csv
@@ -0,0 +1,31 @@
+config,method,parameter,rmse_mean,rmse_sd
+0.26.0.dev1 (Phase 0 commit),cpm-empirical,alpha,"0.013 [0.008, 0.018]","0.026 [0.016, 0.035]"
+0.26.0.dev1 (Phase 0 commit),cpm-empirical,temperature,"0.199 [0.129, 0.262]","0.212 [0.136, 0.281]"
+0.26.0.dev1 (Phase 0 commit),cpm-variational,alpha,"0.012 [0.008, 0.015]","0.028 [0.020, 0.035]"
+0.26.0.dev1 (Phase 0 commit),cpm-variational,temperature,"0.184 [0.122, 0.240]","0.220 [0.134, 0.296]"
+"new, defaults",cpm-empirical,alpha,"0.013 [0.008, 0.018]","0.026 [0.016, 0.035]"
+"new, defaults",cpm-empirical,temperature,"0.199 [0.129, 0.262]","0.212 [0.136, 0.281]"
+"new, defaults",cpm-variational,alpha,"0.012 [0.008, 0.015]","0.028 [0.020, 0.035]"
+"new, defaults",cpm-variational,temperature,"0.184 [0.122, 0.240]","0.220 [0.134, 0.296]"
+"new, hessian",cpm-empirical,alpha,"0.012 [0.008, 0.017]","0.021 [0.011, 0.028]"
+"new, hessian",cpm-empirical,temperature,"0.197 [0.128, 0.261]","0.191 [0.109, 0.256]"
+"new, hessian",cpm-variational,alpha,"0.012 [0.008, 0.015]","0.025 [0.017, 0.031]"
+"new, hessian",cpm-variational,temperature,"0.180 [0.121, 0.235]","0.170 [0.105, 0.225]"
+"new, start",cpm-empirical,alpha,"0.013 [0.008, 0.017]","0.027 [0.016, 0.037]"
+"new, start",cpm-empirical,temperature,"0.199 [0.130, 0.262]","0.211 [0.136, 0.281]"
+"new, start",cpm-variational,alpha,"0.012 [0.009, 0.015]","0.029 [0.021, 0.036]"
+"new, start",cpm-variational,temperature,"0.186 [0.124, 0.242]","0.211 [0.133, 0.282]"
+"new, start + hessian",cpm-empirical,alpha,"0.012 [0.008, 0.017]","0.021 [0.012, 0.028]"
+"new, start + hessian",cpm-empirical,temperature,"0.200 [0.129, 0.262]","0.191 [0.108, 0.256]"
+"new, start + hessian",cpm-variational,alpha,"0.012 [0.008, 0.015]","0.025 [0.017, 0.031]"
+"new, start + hessian",cpm-variational,temperature,"0.183 [0.123, 0.236]","0.170 [0.105, 0.225]"
+"new, start + hessian, 1 start",cpm-empirical,alpha,"0.012 [0.008, 0.016]","0.022 [0.013, 0.028]"
+"new, start + hessian, 1 start",cpm-empirical,temperature,"0.200 [0.130, 0.262]","0.180 [0.094, 0.249]"
+"new, start + hessian, 1 start",cpm-variational,alpha,"0.011 [0.009, 0.014]","0.025 [0.018, 0.031]"
+"new, start + hessian, 1 start",cpm-variational,temperature,"0.183 [0.122, 0.236]","0.171 [0.107, 0.226]"
+comparison (cpm-empirical),cpm-empirical,alpha,"0.013 [0.008, 0.018]","0.026 [0.016, 0.035]"
+comparison (cpm-empirical),cpm-empirical,temperature,"0.199 [0.127, 0.259]","0.213 [0.133, 0.279]"
+comparison (cpm-variational),cpm-variational,alpha,"0.012 [0.008, 0.015]","0.028 [0.021, 0.035]"
+comparison (cpm-variational),cpm-variational,temperature,"0.184 [0.123, 0.240]","0.205 [0.126, 0.279]"
+comparison (cbm),cbm,alpha,"0.023 [0.018, 0.029]","0.022 [0.015, 0.028]"
+comparison (cbm),cbm,temperature,"0.279 [0.183, 0.362]","0.176 [0.115, 0.234]"
diff --git a/benchmarks/results/hierarchical-comparison-runtime.csv b/benchmarks/results/hierarchical-comparison-runtime.csv
new file mode 100644
index 00000000..6e46141b
--- /dev/null
+++ b/benchmarks/results/hierarchical-comparison-runtime.csv
@@ -0,0 +1,16 @@
+config,method,fits,median_seconds,median_ms_per_participant_iteration,median_iterations,stopped_before_50,seconds_vs_cbm,ms_vs_cbm
+0.26.0.dev1 (Phase 0 commit),cpm-empirical,20,90.87562635005452,20.803921870002526,50.0,0.45,9.902121712865382,4.736912106409524
+0.26.0.dev1 (Phase 0 commit),cpm-variational,20,87.99754375009798,20.113323750020935,50.0,0.45,9.588515905135699,4.579667591858573
+"new, defaults",cpm-empirical,20,73.94655685010366,17.220056519983338,50.0,0.45,8.057471905134292,3.9208902394592307
+"new, defaults",cpm-variational,20,72.27422889997251,16.375839510001242,50.0,0.45,7.875249283171419,3.7286677441042935
+"new, hessian",cpm-empirical,20,58.99186499998905,14.11793514001183,42.0,0.6,6.427957096534108,3.2145582116733054
+"new, hessian",cpm-variational,20,52.845357749960385,13.838825339986943,39.0,0.65,5.758212464850953,3.151006800597213
+"new, start",cpm-empirical,20,64.46359169995412,16.2473555700155,41.0,0.55,7.024175311898925,3.6994128211815522
+"new, start",cpm-variational,20,65.32937605003826,15.764868770822067,44.0,0.65,7.118514161114703,3.5895538448518227
+"new, start + hessian",cpm-empirical,20,47.108444100013,13.017018705396913,38.0,0.7,5.1330985646820135,2.9638869959352707
+"new, start + hessian",cpm-variational,20,40.04392045002896,13.229554900014772,25.5,0.85,4.363323699460235,3.012280047958179
+"new, start + hessian, 1 start",cpm-empirical,20,23.012659750063904,6.773201450007036,38.5,0.65,2.5075387860779714,1.5422120957852348
+"new, start + hessian, 1 start",cpm-variational,20,17.833935800008476,6.9956334311127835,21.5,0.85,1.9432471610249897,1.592858351368021
+comparison (cpm-empirical),cpm-empirical,20,81.6190145999426,19.751574295005412,50.0,0.45,8.893489366879983,4.497299691053249
+comparison (cpm-variational),cpm-variational,20,70.75931959995069,14.922597907213236,50.0,0.45,7.710179540323271,3.3977744738449753
+comparison (cbm),cbm,20,9.17738935,4.391874158241759,21.0,0.9,1.0,1.0
diff --git a/benchmarks/results/hierarchical-profile-after.csv b/benchmarks/results/hierarchical-profile-after.csv
new file mode 100644
index 00000000..d97e34c2
--- /dev/null
+++ b/benchmarks/results/hierarchical-profile-after.csv
@@ -0,0 +1,18 @@
+section,item,value
+fit,seconds without profiler,6.1314
+fit,ms per participant and iteration,12.2628
+evaluations per participant and iteration,optimiser,72.198
+evaluations per participant and iteration,hessian,6.992
+evaluations per participant and iteration,evaluate_fit,1.0
+share of profiled time,fmin_l_bfgs_b,0.8817
+share of profiled time,hessian,0.0775
+share of profiled time,evaluate_fit,0.008
+share of profiled time,update_population,0.0007
+share of profiled time,get_lme,0.0001
+share of profiled time,objective (all callers),0.7932
+share of objective,SessionWrapper.run,0.4603
+share of objective,  of which the RLRW model function,0.3393
+share of objective,Wrapper.reset,0.1391
+share of objective,deepcopy,0.0459
+share of objective,loss,0.2029
+share of objective,Parameters.PDF,0.1232
diff --git a/benchmarks/results/hierarchical-profile-before.csv b/benchmarks/results/hierarchical-profile-before.csv
new file mode 100644
index 00000000..a9c573dc
--- /dev/null
+++ b/benchmarks/results/hierarchical-profile-before.csv
@@ -0,0 +1,18 @@
+section,item,value
+fit,seconds without profiler,8.6299
+fit,ms per participant and iteration,17.2597
+evaluations per participant and iteration,optimiser,70.59
+evaluations per participant and iteration,hessian,22.0
+evaluations per participant and iteration,evaluate_fit,2.0
+share of profiled time,fmin_l_bfgs_b,0.7824
+share of profiled time,numerical_hessian,0.1838
+share of profiled time,evaluate_fit,0.0143
+share of profiled time,update_population,0.0006
+share of profiled time,get_lme,0.0002
+share of profiled time,objective (all callers),0.6966
+share of objective,SessionWrapper.run,0.3991
+share of objective,  of which the RLRW model function,0.3183
+share of objective,Wrapper.reset,0.1337
+share of objective,deepcopy,0.0428
+share of objective,loss,0.2729
+share of objective,Parameters.PDF,0.1246
diff --git a/benchmarks/results/hierarchical-profile-cbm.csv b/benchmarks/results/hierarchical-profile-cbm.csv
new file mode 100644
index 00000000..32a5197d
--- /dev/null
+++ b/benchmarks/results/hierarchical-profile-cbm.csv
@@ -0,0 +1,8 @@
+item,value
+participants,100
+iterations,9
+lap_calls_per_participant,85.59
+hbi_calls_per_participant_and_iteration,71.74
+seconds_instrumented,7.21
+lap_seconds_instrumented,2.75
+ms_per_participant_and_iteration_instrumented,4.956
diff --git a/cpm/__version__.py b/cpm/__version__.py
index 4d482beb..1f113989 100644
--- a/cpm/__version__.py
+++ b/cpm/__version__.py
@@ -1 +1 @@
-__version__ = "0.26.0.dev2"
+__version__ = "0.26.0.dev3"
diff --git a/cpm/core/optimisers.py b/cpm/core/optimisers.py
index 64beefea..4ba684bb 100644
--- a/cpm/core/optimisers.py
+++ b/cpm/core/optimisers.py
@@ -12,6 +12,8 @@
     "group_identifier",
     "prepare_data",
     "numerical_hessian",
+    "finite_difference_hessian",
+    "check_hessian",
 ]
 
 
@@ -34,6 +36,71 @@ def numerical_hessian(func=None, params=None, hessian=None):
     return computed_hessian
 
 
+def check_hessian(hessian):
+    """Raise a ValueError unless `hessian` is one of the optimisers' `hessian` options."""
+    if hessian not in ("numdifftools", "finite_differences"):
+        raise ValueError(
+            f'hessian must be "numdifftools" or "finite_differences", not {hessian!r}.'
+        )
+    return hessian
+
+
+def finite_difference_hessian(func, x, lower, upper, step=1e-4):
+    """
+    The Hessian of `func` at `x` by finite differences that stay within the bounds.
+
+    Parameters
+    ----------
+    func : callable
+        The function, of a 1-D array.
+    x : numpy.ndarray
+        The point, within the bounds, normally the optimum.
+    lower, upper : numpy.ndarray
+        The bounds of each element of `x`, possibly infinite.
+    step : float, optional
+        The step, relative to ``max(1, |x|)``. Default is 1e-4.
+
+    Returns
+    -------
+    numpy.ndarray
+        The Hessian matrix.
+
+    Notes
+    -----
+    If every parameter is at least one step from its bounds, the Hessian is
+    computed with second-order differences (Abramowitz and Stegun, 1972,
+    equations 25.3.23 and 25.3.27), which take d² + d + 1 evaluations for d
+    parameters. Otherwise, forward differences step away from the bounds, which
+    take d(d + 1) / 2 + d + 1 evaluations and are first-order accurate.
+    """
+    x = np.asarray(x, dtype=float)
+    d = x.size
+    h = step * np.maximum(1.0, np.abs(x))
+    f0 = func(x)
+    hessian = np.empty((d, d))
+    if np.all((x - h >= lower) & (x + h <= upper)):
+        steps = np.diag(h)
+        plus = [func(x + steps[i]) for i in range(d)]
+        minus = [func(x - steps[i]) for i in range(d)]
+        for i in range(d):
+            hessian[i, i] = (plus[i] - 2 * f0 + minus[i]) / h[i] ** 2
+            for j in range(i + 1, d):
+                both = func(x + steps[i] + steps[j]) + func(x - steps[i] - steps[j])
+                hessian[i, j] = hessian[j, i] = (
+                    both - plus[i] - plus[j] + 2 * f0 - minus[i] - minus[j]
+                ) / (2 * h[i] * h[j])
+        return hessian
+    h = np.where(x + 2 * h > upper, -h, h)
+    steps = np.diag(h)
+    single = [func(x + steps[i]) for i in range(d)]
+    for i in range(d):
+        for j in range(i, d):
+            hessian[i, j] = hessian[j, i] = (
+                func(x + steps[i] + steps[j]) - single[i] - single[j] + f0
+            ) / (h[i] * h[j])
+    return hessian
+
+
 def objective(pars, function, data, loss, prior=False):
     """
     The `objective` function calculates a metric by comparing predicted values with
diff --git a/cpm/hierarchical/_chains.py b/cpm/hierarchical/_chains.py
index 8a9bb2dc..832cc77d 100644
--- a/cpm/hierarchical/_chains.py
+++ b/cpm/hierarchical/_chains.py
@@ -1,8 +1,57 @@
-"""The random starting priors of the chains of the hierarchical methods."""
+"""The random starting priors of the chains of the hierarchical methods, and where their participant-wise fits start."""
 
 import numpy as np
 
 
+def positive_definite(hessians):
+    """
+    Whether each participant's Hessian is positive definite (finite, with a Cholesky decomposition).
+
+    Parameters
+    ----------
+    hessians : array-like
+        The Hessians of the negative log posterior density, shape (participants, parameters, parameters).
+
+    Returns
+    -------
+    numpy.ndarray
+        One boolean per participant.
+    """
+    out = np.zeros(len(hessians), dtype=bool)
+    for i, hessian in enumerate(hessians):
+        if np.all(np.isfinite(hessian)):
+            try:
+                np.linalg.cholesky(hessian)
+                out[i] = True
+            except np.linalg.LinAlgError:
+                pass
+    return out
+
+
+def check_start(start):
+    if start not in ("random", "prior_mean"):
+        raise ValueError(f'start must be "random" or "prior_mean", not {start!r}.')
+    return start
+
+
+def start_from_prior_locations(optimiser):
+    """
+    Replace the first start of `optimiser` with the location of the current prior of each free parameter, clipped into its bounds.
+
+    The other starts are the random guesses that `optimiser.reset()` drew, so the
+    same random numbers are drawn as with random starts only.
+
+    Parameters
+    ----------
+    optimiser : object
+        An optimiser of `cpm.optimisation`, after `reset()`.
+    """
+    parameters = optimiser.model.parameters
+    locations = [getattr(parameters, name).prior.kwds["loc"] for name in parameters.free()]
+    lower, upper = parameters.bounds()
+    optimiser.initial_guess[0] = np.clip(locations, lower, upper)
+
+
 def starting_priors(names, bounds, priors):
     """
     Draw a random starting prior for each parameter, for every chain after the first.
diff --git a/cpm/hierarchical/empirical.py b/cpm/hierarchical/empirical.py
index 2e5f6052..78635a67 100644
--- a/cpm/hierarchical/empirical.py
+++ b/cpm/hierarchical/empirical.py
@@ -2,7 +2,7 @@
 import pandas as pd
 import copy
 from ..core.diagnostics import convergence_diagnostics_plots, parameter_bounds
-from ._chains import starting_priors
+from ._chains import check_start, positive_definite, start_from_prior_locations, starting_priors
 
 
 class EmpiricalBayes:
@@ -23,6 +23,8 @@ class EmpiricalBayes:
         The number of random parameter initialisations. Default is 4. The first chain starts from the priors of the model, and every later chain from priors drawn at random within the bounds of the parameters; `numpy.random.seed` makes them reproducible.
     quiet : bool, optional
         Whether to suppress the output. Default is False.
+    start : {"random", "prior_mean"}, optional
+        Where the participant-wise fits start on every iteration. With "random", every start is drawn at random within the bounds of the parameters. With "prior_mean", the first start is the mean (location) of the current prior of each parameter, clipped into its bounds, as the cbm toolbox does, and the others are drawn at random. Default is "random".
 
     Notes
     -----
@@ -49,6 +51,8 @@ class EmpiricalBayes:
     - When calculating the within-subject variance via the Hessian matrix, the algorithm sets any non-finite or non-positive values to NaN.
     - If the second derivative of the negative log posterior density function is not finite, we set the log determinant to -1e6.
 
+    An estimate on a bound of a parameter can have a Hessian matrix that is not positive definite. The estimate is kept, as the bounds require. The `positive_definite` column of `fit` records whether each participant's Hessian matrix was positive definite on each iteration, and the `not_positive_definite` column of `hyperparameters` how many were not.
+
     References
     ----------
 
@@ -89,8 +93,10 @@ def __init__(
         tolerance=1e-6,
         chain=4,
         quiet=False,
+        start="random",
         **kwargs,
     ):
+        self.start = check_start(start)
         self.function = copy.deepcopy(optimiser.model)
         self.optimiser = copy.deepcopy(optimiser)
         # bounds here should include mean and std for all parameters
@@ -219,6 +225,8 @@ def __log_det_hessian_batch(hessian):
         for iteration in range(self.iteration):
 
             self.optimiser.reset()
+            if self.start == "prior_mean":
+                start_from_prior_locations(self.optimiser)
 
             # perform participant-wise optimisation, extracting MAP parameter estimates,
             # the Hessian matrix of the target function evaluated at the MAP parameter estimates,
@@ -248,6 +256,9 @@ def __log_det_hessian_batch(hessian):
             parameter_long["ppt"] = np.arange(len(parameters))
             parameter_long["iteration"] = iteration + 1
             parameter_long["chain"] = chain_index
+            ## an estimate on a bound can have a Hessian that is not positive definite
+            definite = positive_definite(hessian)
+            parameter_long["positive_definite"] = definite
             fit_records.append(parameter_long)
             # turn any non-finite values into NaN, to avoid subsequent issues
             # with calculating parameter means and variances
@@ -333,6 +344,7 @@ def __log_det_hessian_batch(hessian):
                         "sd": stdev[i],
                         "lme": summed_lme,
                         "reject": reject,
+                        "not_positive_definite": int((~definite).sum()),
                     }
                 )
 
diff --git a/cpm/hierarchical/variational.py b/cpm/hierarchical/variational.py
index dc5a0b28..a0453740 100644
--- a/cpm/hierarchical/variational.py
+++ b/cpm/hierarchical/variational.py
@@ -8,7 +8,7 @@
 
 from ..generators import Parameters
 from ..core.diagnostics import convergence_diagnostics_plots, parameter_bounds
-from ._chains import starting_priors
+from ._chains import check_start, positive_definite, start_from_prior_locations, starting_priors
 
 
 class VariationalBayes:
@@ -33,6 +33,8 @@ class VariationalBayes:
         A dictionary of given parameter values of the prior distributions on the population-level parameters (means mu and precisions tau). See Notes for details. Default is None.
     convergence : str, optional
         The convergence criterion. Default is 'parameters'. Options are 'lme' and 'parameters'.
+    start : {"random", "prior_mean"}, optional
+        Where the participant-wise fits start on every iteration. With "random", every start is drawn at random within the bounds of the parameters. With "prior_mean", the first start is the mean (location) of the current prior of each parameter, clipped into its bounds, as in the HBI implementation of Piray et al. (2019), and the others are drawn at random. Default is "random".
 
     Notes
     -----
@@ -58,6 +60,8 @@ class VariationalBayes:
 
     The convergence criterion can be set to 'lme' or 'parameters'. If set to 'lme', the algorithm will stop when the log model evidence converges. If set to 'parameters', the algorithm will stop when the "normalized" means of the population-level parameters converge.
 
+    An estimate on a bound of a parameter can have a Hessian matrix that is not positive definite. The estimate is kept, as the bounds require. The `positive_definite` column of `fit` records whether each participant's Hessian matrix was positive definite on each iteration, and the `not_positive_definite` column of `hyperparameters` how many were not.
+
     References
     ----------
 
@@ -98,8 +102,10 @@ def __init__(
         hyperpriors=None,
         convergence="parameters",
         quiet=False,
+        start="random",
         **kwargs,
     ):
+        self.start = check_start(start)
         # write input arguments to self
         self.optimiser = copy.deepcopy(optimiser)
         self.objective = objective
@@ -203,6 +209,8 @@ def update_participants(self, iter_idx=0, chain_idx=0):
 
         # delete any pre-existing optimisation output / parameter estimates
         self.optimiser.reset()
+        if self.start == "prior_mean":
+            start_from_prior_locations(self.optimiser)
 
         # run the optimisation
         self.optimiser.optimise()
@@ -226,7 +234,6 @@ def update_participants(self, iter_idx=0, chain_idx=0):
         parameter_long["ppt"] = np.arange(self.__n_ppt__)
         parameter_long["iteration"] = iter_idx + 1
         parameter_long["chain"] = chain_idx
-        self.__fit_chunks__.append(parameter_long)
 
         # extract the participant-wise unnormalised log posterior density at the
         # optimised parameter values
@@ -247,6 +254,10 @@ def update_participants(self, iter_idx=0, chain_idx=0):
         if self.objective != "minimise":
             hessian = -1 * hessian
 
+        ## an estimate on a bound can have a Hessian that is not positive definite
+        parameter_long["positive_definite"] = positive_definite(hessian)
+        self.__fit_chunks__.append(parameter_long)
+
         return param, log_posterior, hessian
 
     def get_lme(self, log_post, hessian):
@@ -499,6 +510,7 @@ def __inv_hessian_batch(hessian):
         # buffered (see the `hyperparameters` property above) rather than
         # concatenated into self.hyperparameters on every single parameter,
         # to avoid an O(n^2) cost across many iterations x parameters.
+        not_definite = int((~positive_definite(hessian)).sum())
         for i, name in enumerate(self.__param_names__):
             self.__hyper_chunks__.append(
                 {
@@ -509,6 +521,7 @@ def __inv_hessian_batch(hessian):
                     "mean_se": E_mu_error[i],
                     "sd": E_sd[i],
                     "lme": lme,
+                    "not_positive_definite": not_definite,
                 }
             )
 
diff --git a/cpm/optimisation/bads.py b/cpm/optimisation/bads.py
index 9b853cc6..5b972fcf 100644
--- a/cpm/optimisation/bads.py
+++ b/cpm/optimisation/bads.py
@@ -5,6 +5,8 @@
     evaluate_fit,
     fit_extras,
     numerical_hessian,
+    finite_difference_hessian,
+    check_hessian,
     group_identifier,
     prepare_data,
 )
@@ -95,6 +97,8 @@ class Bads:
         The libraries required for the parallel processing with `ipyparallel` with the IPython kernel. Default is `["numpy", "pandas"]`.
     ppt_identifier : str
         The column (or the key in the participant data dictionaries) that contains the participant identifier, returned in the optimization details. Required if `data` is a pd.DataFrame. If `data` is a pd.DataFrameGroupBy grouped by a single column, such as `data.groupby("ppt")`, it defaults to the name of that column. Default is `None`.
+    hessian : {"numdifftools", "finite_differences"}
+        How the Hessian matrix recorded in `fit` is computed. With "numdifftools", it is computed with numdifftools' forward differences at the optimised parameter values plus 1e-3, for every start. With "finite_differences", it is computed at the optimised parameter values with finite differences that stay within the bounds of the parameters (see :func:`cpm.core.optimisers.finite_difference_hessian`), once per participant, for the best start only, together with `log_likelihood`, `log_prior` and `metrics`. Default is "numdifftools".
     **kwargs : dict
         Additional keyword arguments. See the `pybads.bads <https://acerbilab.github.io/pybads/api/classes/bads.html>`__ documentation for what is supported.
 
@@ -137,6 +141,7 @@ def __init__(
         libraries=["numpy", "pandas"],
         ppt_identifier=None,
         display=True,
+        hessian="numdifftools",
         **kwargs,
     ):
         self.model = copy.deepcopy(model)
@@ -150,6 +155,7 @@ def __init__(
         self.prior = prior
         self.metrics = metrics
         self.kwargs = numpy_bool_options(kwargs)
+        self.hessian = check_hessian(hessian)
 
         self.fit = []
         self.details = []
@@ -243,6 +249,10 @@ def target(x):
 
             result = optimizer.optimize()
             result = dict(result.items())
+            if self.hessian == "finite_differences":
+                ## the Hessian and the extras follow for the best start (see __complete)
+                result.update({"ppt": ppt})
+                return result
 
             def f(x):
                 return objective(x, model, observed, loss, prior)
@@ -262,16 +272,56 @@ def f(x):
             result.update({"ppt": ppt})
             return result
 
+        def __complete(item):
+            ## the Hessian and the extras of a participant's best start
+            participant, result = item
+            participant_dc, observed, ppt = decompose(
+                participant=participant,
+                pandas=self.__pandas__,
+                identifier=self.ppt_identifier,
+            )
+            model.reset(data=participant_dc)
+
+            def f(x):
+                return objective(x, model, observed, loss, prior)
+
+            result = {**result, "hessian": finite_difference_hessian(f, result["x"], lower, upper)}
+            if prior or self.metrics:
+                result.update(
+                    evaluate_fit(
+                        result["x"], model, observed, loss, prior, self.metrics
+                    )
+                )
+            ## the identifier last, as in __task
+            result["ppt"] = result.pop("ppt")
+            return result
+
         def __extract_nll(result):
             output = np.zeros(len(result))
             for i in range(len(result)):
                 output[i] = result[i].get("fval")
             return output.copy()
 
+        def __record(results):
+            self.details = copy.deepcopy(results)
+            for result in results:
+                parameters = {}
+                for i in range(len(self.parameter_names)):
+                    parameters[self.parameter_names[i]] = copy.deepcopy(
+                        result["x"][i]
+                    )
+
+                self.parameters.append(copy.deepcopy(parameters))
+                self.fit.append(
+                    __unpack(copy.deepcopy(result), id=self.ppt_identifier)
+                )
+
         loss = self.loss
         model = self.model
         prior = self.prior
+        lower, upper = (np.asarray(b, dtype=float) for b in self.__bounds__)
 
+        starts = []
         for i in range(len(self.initial_guess)):
             if self.display:
                 print(
@@ -291,21 +341,11 @@ def __extract_nll(result):
                 results = list(map(__task, self.data))
 
             ## extract the negative log likelihoods for each ppt
-            if i == 0:
-                self.details = copy.deepcopy(results)
+            if self.hessian == "finite_differences":
+                starts.append(results)
+            elif i == 0:
                 old_nll = __extract_nll(results)
-                parameters = {}
-                for result in results:
-                    parameters = {}
-                    for i in range(len(self.parameter_names)):
-                        parameters[self.parameter_names[i]] = copy.deepcopy(
-                            result["x"][i]
-                        )
-
-                    self.parameters.append(copy.deepcopy(parameters))
-                    self.fit.append(
-                        __unpack(copy.deepcopy(result), id=self.ppt_identifier)
-                    )
+                __record(results)
             else:
                 nll = __extract_nll(results)
                 # check if ppt fit is better than the previous fit
@@ -320,6 +360,28 @@ def __extract_nll(result):
                         copy.deepcopy(results[ppt]), id=self.ppt_identifier
                     )
 
+        if self.hessian == "finite_differences":
+            ## choose each participant's start as above: a later start replaces
+            ## the first where it is better
+            best = list(starts[0])
+            first_nll = __extract_nll(starts[0])
+            for results in starts[1:]:
+                for ppt in np.where(__extract_nll(results) < first_nll)[0]:
+                    best[ppt] = results[ppt]
+            items = list(zip(self.data, best))
+            if self.__parallel__:
+                results = execute_parallel(
+                    job=__complete,
+                    data=items,
+                    method=None,
+                    pandas=False,
+                    cl=self.cl,
+                    libraries=self.__libraries__,
+                )
+            else:
+                results = list(map(__complete, items))
+            __record(results)
+
         return None
 
     def reset(self, initial_guess=True):
diff --git a/cpm/optimisation/fmin.py b/cpm/optimisation/fmin.py
index 6f445d9d..a916cf83 100644
--- a/cpm/optimisation/fmin.py
+++ b/cpm/optimisation/fmin.py
@@ -4,6 +4,8 @@
     evaluate_fit,
     fit_extras,
     numerical_hessian,
+    finite_difference_hessian,
+    check_hessian,
     group_identifier,
     prepare_data,
 )
@@ -27,6 +29,61 @@
 LBFGSB_VERBOSITY_SUPPORTED = "disp" in inspect.signature(fmin_l_bfgs_b).parameters
 
 
+def gradient_by_forward_differences(fun, x, f0, lower, upper, epsilon):
+    """
+    The gradient of `fun` at `x` by forward differences, as SciPy's L-BFGS-B estimates it with `approx_grad=True`.
+
+    Parameters
+    ----------
+    fun : callable
+        The function, of a 1-D array.
+    x : numpy.ndarray
+        The point, within the bounds.
+    f0 : float
+        `fun(x)`.
+    lower, upper : numpy.ndarray
+        The bounds of each element of `x`, possibly infinite.
+    epsilon : float
+        The step.
+
+    Returns
+    -------
+    numpy.ndarray
+        The gradient.
+
+    Notes
+    -----
+    This follows `scipy.optimize._numdiff.approx_derivative` with the "2-point"
+    method and an absolute step: a step that would leave the bounds is taken in
+    the other direction, and each difference is divided by the step as it is
+    represented, ``(x + h) - x``. The gradient is therefore exactly the one SciPy
+    computes, so the optimiser takes exactly the same path, but without SciPy's
+    overhead on every evaluation. It works on Python floats, which follow the
+    same IEEE arithmetic as NumPy's and are faster for a few parameters.
+    """
+    gradient = np.empty(x.size)
+    for i, (value, low, high) in enumerate(zip(x.tolist(), lower.tolist(), upper.tolist())):
+        h = epsilon
+        ## a step too small to change x falls back to SciPy's default relative step
+        if (value + h) - value == 0:
+            h = _SQRT_EPS * (1.0 if value >= 0 else -1.0) * max(1.0, abs(value))
+        lower_distance, upper_distance = value - low, high - value
+        if abs(h) <= max(lower_distance, upper_distance):
+            if value + h < low or value + h > high:
+                h = -h
+        elif upper_distance >= lower_distance:
+            h = upper_distance
+        else:
+            h = -lower_distance
+        step = x.copy()
+        step[i] = value + h
+        gradient[i] = (fun(step) - f0) / ((value + h) - value)
+    return gradient
+
+
+_SQRT_EPS = float(np.finfo(float).eps ** 0.5)
+
+
 def lbfgsb_options(display=False, **kwargs):
     """
     Assemble the keyword arguments for `scipy.optimize.fmin_l_bfgs_b`, keeping the solver's verbosity options only where the installed SciPy still accepts them.
@@ -98,6 +155,8 @@ class Fmin:
         The libraries to import for parallel processing for `ipyparallel` with the IPython kernel. Default is `["numpy", "pandas"]`.
     ppt_identifier : str
         The column (or the key in the participant data dictionaries) that contains the participant identifier, returned in the optimization details. Required if `data` is a pd.DataFrame. If `data` is a pd.DataFrameGroupBy grouped by a single column, such as `data.groupby("ppt")`, it defaults to the name of that column. Default is `None`.
+    hessian : {"numdifftools", "finite_differences"}
+        How the Hessian matrix recorded in `fit` is computed. With "numdifftools", it is computed with numdifftools' forward differences at the optimised parameter values plus 1e-3, for every start. With "finite_differences", it is computed at the optimised parameter values with finite differences that stay within the bounds of the parameters, if the optimum is within them (see :func:`cpm.core.optimisers.finite_difference_hessian`), once per participant, for the best start only, together with `log_likelihood`, `log_prior` and `metrics`. Default is "numdifftools".
     **kwargs : dict
         Additional keyword arguments. See the :func:`scipy.optimize.fmin` documentation for what is supported.
 
@@ -136,6 +195,7 @@ def __init__(
         number_of_starts=1,
         ppt_identifier=None,
         display=False,
+        hessian="numdifftools",
         **kwargs,
     ):
         self.model = copy.deepcopy(model)
@@ -150,6 +210,7 @@ def __init__(
         self.metrics = metrics
         self.kwargs = kwargs
         self.display = display
+        self.hessian = check_hessian(hessian)
 
         self.fit = []
         self.details = []
@@ -222,6 +283,9 @@ def __task(participant, **args):
                 **self.kwargs,
                 full_output=True,
             )
+            if self.hessian == "finite_differences":
+                ## the Hessian and the extras follow for the best start (see __complete)
+                return (*result, ppt)
 
             def f(x):
                 return objective(x, model, observed, loss, prior)
@@ -241,16 +305,54 @@ def f(x):
             result = (*result, ppt)
             return result
 
+        def __complete(item):
+            ## the Hessian and the extras of a participant's best start
+            participant, result = item
+            participant_dc, observed, ppt = decompose(
+                participant=participant,
+                pandas=self.__pandas__,
+                identifier=self.ppt_identifier,
+            )
+            model.reset(data=participant_dc)
+
+            def f(x):
+                return objective(x, model, observed, loss, prior)
+
+            result = (*result[:-1], finite_difference_hessian(f, result[0], lower, upper))
+            if prior or self.metrics:
+                result = (
+                    *result,
+                    *evaluate_fit(
+                        result[0], model, observed, loss, prior, self.metrics
+                    ).values(),
+                )
+            return (*result, ppt)
+
         def __extract_nll(result):
             output = np.zeros(len(result))
             for i in range(len(result)):
                 output[i] = result[i][1]
             return output.copy()
 
+        def __record(results):
+            self.details = copy.deepcopy(results)
+            parameters = {}
+            for result in results:
+                for i in range(len(self.parameter_names)):
+                    parameters[self.parameter_names[i]] = copy.deepcopy(
+                        result[0][i]
+                    )
+                self.parameters.append(copy.deepcopy(parameters))
+                self.fit.append(
+                    __unpack(copy.deepcopy(result), id=self.ppt_identifier)
+                )
+
         loss = self.loss
         model = self.model
         prior = self.prior
+        lower, upper = (np.asarray(b, dtype=float) for b in self.model.parameters.bounds())
 
+        starts = []
         for i in range(len(self.initial_guess)):
             if self.display:
                 print(
@@ -270,19 +372,11 @@ def __extract_nll(result):
                 results = list(map(__task, self.data))
 
             ## extract the negative log likelihoods for each ppt
-            if i == 0:
+            if self.hessian == "finite_differences":
+                starts.append(results)
+            elif i == 0:
                 old_nll = __extract_nll(results)
-                self.details = copy.deepcopy(results)
-                parameters = {}
-                for result in results:
-                    for i in range(len(self.parameter_names)):
-                        parameters[self.parameter_names[i]] = copy.deepcopy(
-                            result[0][i]
-                        )
-                    self.parameters.append(copy.deepcopy(parameters))
-                    self.fit.append(
-                        __unpack(copy.deepcopy(result), id=self.ppt_identifier)
-                    )
+                __record(results)
             else:
                 nll = __extract_nll(results)
                 # check if ppt fit is better than the previous fit
@@ -297,6 +391,28 @@ def __extract_nll(result):
                         copy.deepcopy(results[ppt]), id=self.ppt_identifier
                     )
 
+        if self.hessian == "finite_differences":
+            ## choose each participant's start as above: a later start replaces
+            ## the first where it is better
+            best = list(starts[0])
+            first_nll = __extract_nll(starts[0])
+            for results in starts[1:]:
+                for ppt in np.where(__extract_nll(results) < first_nll)[0]:
+                    best[ppt] = results[ppt]
+            items = list(zip(self.data, best))
+            if self.__parallel__:
+                results = execute_parallel(
+                    job=__complete,
+                    data=items,
+                    pandas=False,
+                    method=None,
+                    cl=self.cl,
+                    libraries=self.__libraries__,
+                )
+            else:
+                results = list(map(__complete, items))
+            __record(results)
+
         return None
 
     def reset(self, initial_guess=True):
@@ -370,6 +486,8 @@ class FminBound:
         The libraries to import for parallel processing for `ipyparallel` with the IPython kernel. Default is `["numpy", "pandas"]`.
     ppt_identifier : str
         The column (or the key in the participant data dictionaries) that contains the participant identifier, returned in the optimization details. Required if `data` is a pd.DataFrame. If `data` is a pd.DataFrameGroupBy grouped by a single column, such as `data.groupby("ppt")`, it defaults to the name of that column. Default is `None`.
+    hessian : {"numdifftools", "finite_differences"}
+        How the Hessian matrix recorded in `fit` is computed. With "numdifftools", it is computed with numdifftools' forward differences at the optimised parameter values plus 1e-3, for every start; at an upper bound, that point is outside the bounds, which can give a Hessian of zeros. With "finite_differences", it is computed at the optimised parameter values with finite differences that stay within the bounds (see :func:`cpm.core.optimisers.finite_difference_hessian`), once per participant, for the best start only, together with `log_likelihood`, `log_prior` and `metrics`. Default is "numdifftools".
     **kwargs : dict
         Additional keyword arguments. See the :func:`scipy.optimize.fmin_l_bfgs_b` documentation for what is supported. The solver's own `disp` and `iprint` options, which `display` also sets, were removed in SciPy 1.18.0 and are ignored there.
 
@@ -408,6 +526,7 @@ def __init__(
         metrics=None,
         ppt_identifier=None,
         display=False,
+        hessian="numdifftools",
         **kwargs,
     ):
         self.model = copy.deepcopy(model)
@@ -422,6 +541,7 @@ def __init__(
         self.metrics = metrics
         self.kwargs = kwargs
         self.display = display
+        self.hessian = check_hessian(hessian)
 
         self.fit = []
         self.details = []
@@ -483,6 +603,17 @@ def __unpack(x, id=None):
         model = self.model
         prior = self.prior
         options = lbfgsb_options(display=self.display, **self.kwargs)
+        lower, upper = np.asarray(bounds, dtype=float).T
+        ## with `approx_grad`, the gradient is computed here, exactly as SciPy
+        ## would, together with the objective (see `gradient_by_forward_differences`)
+        approx_grad = options.pop("approx_grad", False)
+        if approx_grad:
+            options.pop("fprime", None)
+            epsilon = options.pop("epsilon", 1e-8)
+            evaluations_per_point = len(self.parameter_names) + 1
+            ## SciPy counts the evaluations of the gradient towards `maxfun`, but
+            ## here it only counts the points, so the limit is in points
+            options["maxfun"] = options.get("maxfun", 15000) // evaluations_per_point
 
         def __task(participant, **args):
 
@@ -494,20 +625,39 @@ def __task(participant, **args):
 
             model.reset(data=participant_dc)
 
-            result = fmin_l_bfgs_b(
-                objective,
-                x0=self.__current_guess__,
-                bounds=bounds,
-                args=(model, observed, loss, prior),
-                **options,
-            )
-
             def f(x):
                 return objective(x, model, observed, loss, prior)
 
+            if approx_grad:
+
+                def objective_and_gradient(x):
+                    f0 = f(x)
+                    return f0, gradient_by_forward_differences(f, x, f0, lower, upper, epsilon)
+
+                result = fmin_l_bfgs_b(
+                    objective_and_gradient,
+                    x0=self.__current_guess__,
+                    bounds=bounds,
+                    **options,
+                )
+                result[2]["funcalls"] *= evaluations_per_point
+            else:
+                result = fmin_l_bfgs_b(
+                    objective,
+                    x0=self.__current_guess__,
+                    bounds=bounds,
+                    args=(model, observed, loss, prior),
+                    **options,
+                )
+
+            result = (*result[0:2], *tuple(list(result[2].values())))
+            if self.hessian == "finite_differences":
+                ## the Hessian and the extras follow for the best start (see __complete)
+                return (*result, ppt)
+
             hessian = numerical_hessian(func=f, params=result[0] + 1e-3)
 
-            result = (*result[0:2], *tuple(list(result[2].values())), hessian)
+            result = (*result, hessian)
 
             if prior or self.metrics:
                 result = (
@@ -520,12 +670,49 @@ def f(x):
             result = (*result, ppt)
             return result
 
+        def __complete(item):
+            ## the Hessian and the extras of a participant's best start
+            participant, result = item
+            participant_dc, observed, ppt = decompose(
+                participant=participant,
+                pandas=self.__pandas__,
+                identifier=self.ppt_identifier,
+            )
+            model.reset(data=participant_dc)
+
+            def f(x):
+                return objective(x, model, observed, loss, prior)
+
+            result = (*result[:-1], finite_difference_hessian(f, result[0], lower, upper))
+            if prior or self.metrics:
+                result = (
+                    *result,
+                    *evaluate_fit(
+                        result[0], model, observed, loss, prior, self.metrics
+                    ).values(),
+                )
+            return (*result, ppt)
+
         def __extract_nll(result):
             output = np.zeros(len(result))
             for i in range(len(result)):
                 output[i] = result[i][1]
             return output.copy()
 
+        def __record(results):
+            self.details = copy.deepcopy(results)
+            parameters = {}
+            for result in results:
+                for i in range(len(self.parameter_names)):
+                    parameters[self.parameter_names[i]] = copy.deepcopy(
+                        result[0][i]
+                    )
+                self.parameters.append(copy.deepcopy(parameters))
+                self.fit.append(
+                    __unpack(copy.deepcopy(result), id=self.ppt_identifier)
+                )
+
+        starts = []
         for i in range(len(self.initial_guess)):
             if self.display:
                 print(
@@ -545,19 +732,11 @@ def __extract_nll(result):
                 results = list(map(__task, self.data))
 
             ## extract the negative log likelihoods for each ppt
-            if i == 0:
+            if self.hessian == "finite_differences":
+                starts.append(results)
+            elif i == 0:
                 old_nll = __extract_nll(results)
-                self.details = copy.deepcopy(results)
-                parameters = {}
-                for result in results:
-                    for i in range(len(self.parameter_names)):
-                        parameters[self.parameter_names[i]] = copy.deepcopy(
-                            result[0][i]
-                        )
-                    self.parameters.append(copy.deepcopy(parameters))
-                    self.fit.append(
-                        __unpack(copy.deepcopy(result), id=self.ppt_identifier)
-                    )
+                __record(results)
             else:
                 nll = __extract_nll(results)
                 # check if ppt fit is better than the previous fit
@@ -572,6 +751,28 @@ def __extract_nll(result):
                         copy.deepcopy(results[ppt]), id=self.ppt_identifier
                     )
 
+        if self.hessian == "finite_differences":
+            ## choose each participant's start as above: a later start replaces
+            ## the first where it is better
+            best = list(starts[0])
+            first_nll = __extract_nll(starts[0])
+            for results in starts[1:]:
+                for ppt in np.where(__extract_nll(results) < first_nll)[0]:
+                    best[ppt] = results[ppt]
+            items = list(zip(self.data, best))
+            if self.__parallel__:
+                results = execute_parallel(
+                    job=__complete,
+                    data=items,
+                    method=None,
+                    cl=self.cl,
+                    pandas=False,
+                    libraries=self.__libraries__,
+                )
+            else:
+                results = list(map(__complete, items))
+            __record(results)
+
         return None
 
     def reset(self, initial_guess=True):
diff --git a/cpm/optimisation/free.py b/cpm/optimisation/free.py
index 39eaf514..304492ad 100644
--- a/cpm/optimisation/free.py
+++ b/cpm/optimisation/free.py
@@ -4,6 +4,8 @@
     evaluate_fit,
     fit_extras,
     numerical_hessian,
+    finite_difference_hessian,
+    check_hessian,
     group_identifier,
     prepare_data,
 )
@@ -48,6 +50,8 @@ class Minimize:
         The libraries to import for parallel processing for `ipyparallel` with the IPython kernel. Default is `["numpy", "pandas"]`
     ppt_identifier : str
         The column (or the key in the participant data dictionaries) that contains the participant identifier, returned in the optimization details. Required if `data` is a pd.DataFrame. If `data` is a pd.DataFrameGroupBy grouped by a single column, such as `data.groupby("ppt")`, it defaults to the name of that column. Default is `None`.
+    hessian : {"numdifftools", "finite_differences"}
+        How the Hessian matrix recorded in `details` is computed. With "numdifftools", it is computed with numdifftools' forward differences at the optimised parameter values plus 1e-3, for every start. With "finite_differences", it is computed at the optimised parameter values with finite differences that stay within the bounds of the parameters (see :func:`cpm.core.optimisers.finite_difference_hessian`), once per participant, for the best start only, together with `log_likelihood`, `log_prior` and `metrics`. Default is "numdifftools".
     **kwargs : dict
         Additional keyword arguments. See the :func:`scipy.optimize.minimize` documentation for what is supported.
 
@@ -87,6 +91,7 @@ def __init__(
         number_of_starts=1,
         ppt_identifier=None,
         display=False,
+        hessian="numdifftools",
         **kwargs,
     ):
         self.model = copy.deepcopy(model)
@@ -101,6 +106,7 @@ def __init__(
         self.metrics = metrics
         self.kwargs = kwargs
         self.display = display
+        self.hessian = check_hessian(hessian)
 
         self.method = method
         if isinstance(model, Wrapper):
@@ -173,6 +179,11 @@ def __task(participant, **args):
                 **self.kwargs,
             )
             result = {**result}
+            if self.hessian == "finite_differences":
+                ## the Hessian and the extras follow for the best start (see __complete)
+                if self.ppt_identifier is not None:
+                    result.update({self.ppt_identifier: ppt})
+                return result
 
             def f(x):
                 return objective(x, model, observed, loss, prior)
@@ -193,12 +204,48 @@ def f(x):
                 result.update({self.ppt_identifier: ppt})
             return result
 
+        def __complete(item):
+            ## the Hessian and the extras of a participant's best start
+            participant, result = item
+            participant_dc, observed, ppt = decompose(
+                participant=participant,
+                pandas=self.__pandas__,
+                identifier=self.ppt_identifier,
+            )
+            model.reset(data=participant_dc)
+
+            def f(x):
+                return objective(x, model, observed, loss, prior)
+
+            result = {**result, "hessian": finite_difference_hessian(f, result["x"], lower, upper)}
+            if prior or self.metrics:
+                result.update(
+                    evaluate_fit(
+                        result["x"], model, observed, loss, prior, self.metrics
+                    )
+                )
+            return result
+
         def __extract_nll(result):
             output = np.zeros(len(result))
             for i in range(len(result)):
                 output[i] = result[i]["fun"]
             return output.copy()
 
+        def __record(results):
+            self.details = copy.deepcopy(results)
+            parameters = {}
+            for result in results:
+                for i in range(len(self.parameter_names)):
+                    parameters[self.parameter_names[i]] = copy.deepcopy(
+                        result["x"][i]
+                    )
+                self.parameters.append(copy.deepcopy(parameters))
+                self.fit.append(
+                    __unpack(copy.deepcopy(result), id=self.ppt_identifier)
+                )
+
+        lower, upper = (np.asarray(b, dtype=float) for b in self.model.parameters.bounds())
         bounds = self.model.parameters.bounds()
         bounds = np.asarray(bounds).T
         bounds = list(map(tuple, bounds))
@@ -208,6 +255,7 @@ def __extract_nll(result):
         model = self.model
         prior = self.prior
 
+        starts = []
         for i in range(len(self.initial_guess)):
             if self.display:
                 print(
@@ -227,19 +275,11 @@ def __extract_nll(result):
                 results = list(map(__task, self.data))
 
             ## extract the negative log likelihoods for each ppt
-            if i == 0:
+            if self.hessian == "finite_differences":
+                starts.append(results)
+            elif i == 0:
                 old_nll = __extract_nll(results)
-                self.details = copy.deepcopy(results)
-                parameters = {}
-                for result in results:
-                    for i in range(len(self.parameter_names)):
-                        parameters[self.parameter_names[i]] = copy.deepcopy(
-                            result["x"][i]
-                        )
-                    self.parameters.append(copy.deepcopy(parameters))
-                    self.fit.append(
-                        __unpack(copy.deepcopy(result), id=self.ppt_identifier)
-                    )
+                __record(results)
             else:
                 nll = __extract_nll(results)
                 # check if ppt fit is better than the previous fit
@@ -254,6 +294,28 @@ def __extract_nll(result):
                         copy.deepcopy(results[ppt]), id=self.ppt_identifier
                     )
 
+        if self.hessian == "finite_differences":
+            ## choose each participant's start as above: a later start replaces
+            ## the first where it is better
+            best = list(starts[0])
+            first_nll = __extract_nll(starts[0])
+            for results in starts[1:]:
+                for ppt in np.where(__extract_nll(results) < first_nll)[0]:
+                    best[ppt] = results[ppt]
+            items = list(zip(self.data, best))
+            if self.__parallel__:
+                results = execute_parallel(
+                    job=__complete,
+                    data=items,
+                    method=None,
+                    cl=self.cl,
+                    pandas=False,
+                    libraries=self.__libraries__,
+                )
+            else:
+                results = list(map(__complete, items))
+            __record(results)
+
         return None
 
     def reset(self, initial_guess=True):
diff --git a/cpm/optimisation/minimise.py b/cpm/optimisation/minimise.py
index a979f7f0..ab4babbe 100644
--- a/cpm/optimisation/minimise.py
+++ b/cpm/optimisation/minimise.py
@@ -198,7 +198,9 @@ def bernoulli(predicted=None, observed=None, negative=True, **kwargs):
         np.clip(probabilities, 1e-10, 1 - 1e-10, out=probabilities)
 
         LL = _bernoulli_logpmf(observed.flatten(), probabilities)
-        LL = check_nan_bounds_in_log(LL, bound=bound)
+        ## LL is a new array, which only needs replacing where it is not finite
+        if not np.isfinite(LL).all():
+            LL = check_nan_bounds_in_log(LL, bound=bound)
         LL = np.sum(LL)
         if negative:
             LL = -1 * LL
diff --git a/test/hierarchical/test_positive_definite.py b/test/hierarchical/test_positive_definite.py
new file mode 100644
index 00000000..76c3dc81
--- /dev/null
+++ b/test/hierarchical/test_positive_definite.py
@@ -0,0 +1,73 @@
+"""
+The hierarchical methods record, per participant and iteration, whether the
+Hessian of the negative log posterior is positive definite, and per iteration how
+many are not. An estimate on a bound can have a Hessian that is not; it is kept.
+"""
+
+import warnings
+
+import numpy as np
+import pytest
+
+from cpm.hierarchical import EmpiricalBayes, VariationalBayes
+from cpm.optimisation import FminBound, minimise
+from test.optimisation.test_fminbound_gradient import bandit, rlrw
+
+METHODS = [EmpiricalBayes, VariationalBayes]
+
+
+def method(cls, hessian):
+    warnings.simplefilter("ignore")
+    data = bandit(participants=6, trials=40)
+    ## with an upper bound of 1 on the temperature, estimates end on it, where the
+    ## "numdifftools" Hessian is all zeros (with 0.5, EmpiricalBayes fails on
+    ## these all-zero Hessians, as before)
+    optimiser = FminBound(model=rlrw(data, temperature_upper=1.0), data=data.groupby("ppt"),
+                          minimisation=minimise.LogLikelihood.bernoulli, prior=True, number_of_starts=1,
+                          approx_grad=True, hessian=hessian)
+    tolerance = {"tolerance": 0.0} if cls is EmpiricalBayes else {"tolerance_lme": 0.0, "tolerance_param": 0.0}
+    fit = cls(optimiser=optimiser, iteration=3, chain=1, quiet=True, **tolerance)
+    hessians = []
+    optimise = fit.optimiser.optimise
+
+    def recorded(*args, **kwargs):
+        result = optimise(*args, **kwargs)
+        hessians.append([r["hessian"] for r in fit.optimiser.fit])
+        return result
+
+    fit.optimiser.optimise = recorded
+    np.random.seed(2)
+    fit.optimise()
+    return fit, hessians
+
+
+def is_positive_definite(h):
+    if not np.all(np.isfinite(h)):
+        return False
+    try:
+        np.linalg.cholesky(h)
+        return True
+    except np.linalg.LinAlgError:
+        return False
+
+
+@pytest.mark.parametrize("hessian", ["numdifftools", "finite_differences"])
+@pytest.mark.parametrize("cls", METHODS)
+def test_flags_and_counts(cls, hessian):
+    fit, hessians = method(cls, hessian)
+    expected = np.array([[is_positive_definite(h) for h in step] for step in hessians])
+    flags = fit.fit.sort_values(["iteration", "ppt"]).positive_definite.to_numpy().reshape(expected.shape)
+    np.testing.assert_array_equal(flags, expected)
+    counts = fit.hyperparameters.groupby("iteration").not_positive_definite.first().to_numpy()
+    np.testing.assert_array_equal(counts, (~expected).sum(axis=1))
+    if hessian == "numdifftools":
+        assert counts.sum() > 0  # the all-zero Hessians on the bound
+
+
+@pytest.mark.parametrize("cls", METHODS)
+def test_estimates_on_a_bound_are_kept(cls):
+    fit, _ = method(cls, "numdifftools")
+    flagged = fit.fit[~fit.fit.positive_definite]
+    assert len(flagged) > 0
+    ## each with a parameter on a bound: alpha in [0.01, 1], temperature in [0, 1]
+    assert (flagged.alpha.isin([0.01, 1.0]) | flagged.temperature.isin([0.0, 1.0])).all()
diff --git a/test/hierarchical/test_start.py b/test/hierarchical/test_start.py
new file mode 100644
index 00000000..729a885b
--- /dev/null
+++ b/test/hierarchical/test_start.py
@@ -0,0 +1,91 @@
+"""`start="prior_mean"`: the first start of every participant is the location of the current prior."""
+
+import warnings
+
+import numpy as np
+import pytest
+
+from cpm.applications.reinforcement_learning import RLRW
+from cpm.datasets import load_bandit_data
+from cpm.hierarchical import EmpiricalBayes, VariationalBayes
+from cpm.optimisation import FminBound, minimise
+
+METHODS = [EmpiricalBayes, VariationalBayes]
+
+
+def method(cls, participants=4, iterations=2, chain=1, **kwargs):
+    warnings.simplefilter("ignore")
+    data = load_bandit_data()
+    data["observed"] = data["response"]
+    data = data[data.ppt.isin(sorted(data.ppt.unique())[:participants])].copy()
+    model = RLRW(data=data[data.ppt == data.ppt.iloc[0]], dimensions=4, parameters_settings=[[0.5, 1e-2, 1.0], [2.0, 0.0, 10.0]])
+    model.parameters.update_prior(alpha={"mean": 0.3, "sd": 0.5}, temperature={"mean": 12.0, "sd": 5.0})
+    optimiser = FminBound(
+        model=model,
+        data=data,
+        minimisation=minimise.LogLikelihood.bernoulli,
+        prior=True,
+        ppt_identifier="ppt",
+        number_of_starts=2,
+        approx_grad=True,
+    )
+    tolerance = {"tolerance": 0.0} if cls is EmpiricalBayes else {"tolerance_lme": 0.0, "tolerance_param": 0.0}
+    return cls(optimiser=optimiser, iteration=iterations, chain=chain, quiet=True, **tolerance, **kwargs)
+
+
+def record_guesses(fit):
+    """The guesses and the prior locations the optimiser starts from, on every participant-wise step."""
+    guesses, locations = [], []
+    optimise = fit.optimiser.optimise
+
+    def recorded(*args, **kwargs):
+        guesses.append(fit.optimiser.initial_guess.copy())
+        parameters = fit.optimiser.model.parameters
+        locations.append([getattr(parameters, name).prior.kwds["loc"] for name in parameters.free()])
+        return optimise(*args, **kwargs)
+
+    fit.optimiser.optimise = recorded
+    return guesses, locations
+
+
+@pytest.mark.parametrize("cls", METHODS)
+def test_first_start_is_the_prior_location_clipped_into_the_bounds(cls):
+    fit = method(cls, iterations=3, chain=2, start="prior_mean")
+    guesses, locations = record_guesses(fit)
+    np.random.seed(1)
+    fit.optimise()
+    assert len(guesses) == 6
+    lower, upper = fit.optimiser.model.parameters.bounds()
+    np.testing.assert_array_equal(guesses[0][0], [0.3, 10.0])  # 12 is above the upper bound of 10
+    for guess, location in zip(guesses, locations):
+        np.testing.assert_array_equal(guess[0], np.clip(location, lower, upper))
+
+
+@pytest.mark.parametrize("cls", METHODS)
+def test_other_starts_are_the_random_ones(cls):
+    random, prior_mean = method(cls, start="random"), method(cls, start="prior_mean")
+    random_guesses, _ = record_guesses(random)
+    prior_guesses, _ = record_guesses(prior_mean)
+    np.random.seed(1)
+    random.optimise()
+    np.random.seed(1)
+    prior_mean.optimise()
+    np.testing.assert_array_equal(random_guesses[0][1:], prior_guesses[0][1:])
+    assert not np.array_equal(random_guesses[0][0], prior_guesses[0][0])
+
+
+@pytest.mark.parametrize("cls", METHODS)
+def test_random_is_the_default(cls):
+    default, random = method(cls), method(cls, start="random")
+    np.random.seed(1)
+    default.optimise()
+    np.random.seed(1)
+    random.optimise()
+    assert default.hyperparameters.equals(random.hyperparameters)
+    assert default.fit.equals(random.fit)
+
+
+@pytest.mark.parametrize("cls", METHODS)
+def test_unknown_start_is_an_error(cls):
+    with pytest.raises(ValueError, match="start"):
+        method(cls, start="previous")
diff --git a/test/optimisation/test_fminbound_gradient.py b/test/optimisation/test_fminbound_gradient.py
new file mode 100644
index 00000000..845cfac0
--- /dev/null
+++ b/test/optimisation/test_fminbound_gradient.py
@@ -0,0 +1,131 @@
+"""
+`FminBound` with `approx_grad=True` computes the finite-difference gradient
+itself, together with the objective, instead of leaving it to SciPy. These tests
+check that the fits are exactly the ones SciPy's own `approx_grad` gives.
+"""
+
+import numpy as np
+import pandas as pd
+import pytest
+from scipy.optimize import fmin_l_bfgs_b
+
+from cpm.applications.reinforcement_learning import RLRW
+from cpm.core.data import decompose
+from cpm.core.optimisers import objective
+from cpm.optimisation import FminBound, minimise
+
+LOSS = minimise.LogLikelihood.bernoulli
+
+
+def bandit(participants=4, trials=40, seed=1):
+    rng = np.random.default_rng(seed)
+    arms = np.array([rng.choice(4, size=2, replace=False) for _ in range(participants * trials)]) + 1
+    return pd.DataFrame(
+        {
+            "ppt": np.repeat(np.arange(1, participants + 1), trials),
+            "arm_left": arms[:, 0],
+            "arm_right": arms[:, 1],
+            "reward_left": rng.integers(0, 2, participants * trials).astype(float),
+            "reward_right": rng.integers(0, 2, participants * trials).astype(float),
+            "response": rng.integers(0, 2, participants * trials),
+            "observed": rng.integers(0, 2, participants * trials),
+        }
+    )
+
+
+def rlrw(data, temperature_upper=10.0):
+    model = RLRW(
+        data=data[data.ppt == 1],
+        dimensions=4,
+        parameters_settings=[[0.5, 1e-2, 1.0], [2.0, 0.0, temperature_upper]],
+    )
+    model.parameters.update_prior(alpha={"mean": 0.5, "sd": 0.5}, temperature={"mean": 5.0, "sd": 5.0})
+    return model
+
+
+def scipy_reference(model, data, guess, prior, **options):
+    """What `fmin_l_bfgs_b` with SciPy's own `approx_grad` gives for every participant."""
+    bounds = list(map(tuple, np.asarray(model.parameters.bounds()).T))
+    out = []
+    for participant in data.groupby("ppt"):
+        participant_dc, observed, _ = decompose(participant=participant, pandas=True, identifier="ppt")
+        model.reset(data=participant_dc)
+        out.append(
+            fmin_l_bfgs_b(objective, x0=guess, bounds=bounds, args=(model, observed, LOSS, prior), approx_grad=True, **options)
+        )
+    return out
+
+
+@pytest.mark.parametrize("prior", [False, True])
+@pytest.mark.parametrize("temperature_upper", [10.0, 0.5])
+@pytest.mark.parametrize("options", [{}, {"epsilon": 1e-6}, {"maxfun": 10}])
+def test_fits_are_identical_to_scipy_approx_grad(prior, temperature_upper, options):
+    data = bandit()
+    guess = np.array([0.3, 0.4])
+    fit = FminBound(
+        model=rlrw(data, temperature_upper),
+        data=data.groupby("ppt"),
+        minimisation=LOSS,
+        prior=prior,
+        initial_guess=guess,
+        number_of_starts=1,
+        approx_grad=True,
+        **options,
+    )
+    fit.optimise()
+    reference = scipy_reference(rlrw(data, temperature_upper), data, guess, prior, **options)
+    for got, (x, f, info) in zip(fit.fit, reference):
+        np.testing.assert_array_equal(got["x"], x)
+        assert got["fun"] == f
+        np.testing.assert_array_equal(got["grad"], info["grad"])
+        assert got["funcalls"] == info["funcalls"]
+        assert got["nit"] == info["nit"]
+        assert got["warnflag"] == info["warnflag"]
+        assert got["task"] == info["task"]
+
+
+def test_a_participant_on_the_upper_bound_is_covered():
+    """With an upper bound of 0.5 on the temperature, the estimates end on it, so the step must flip sign."""
+    data = bandit()
+    reference = scipy_reference(rlrw(data, 0.5), data, np.array([0.3, 0.4]), True)
+    assert any(x[1] == 0.5 for x, _, _ in reference)
+
+
+@pytest.mark.parametrize(
+    "x, lower, upper",
+    [
+        ([0.4, 3.0], [1e-2, 0.0], [1.0, 10.0]),  # inside
+        ([1e-2, 10.0], [1e-2, 0.0], [1.0, 10.0]),  # on a lower and an upper bound
+        ([0.4, -3.0], [-np.inf, -np.inf], [np.inf, np.inf]),  # unbounded
+        ([0.5, 1.0], [0.5, 1.0 - 4e-9], [0.5 + 3e-9, 1.0]),  # bounds closer than the step
+        ([1e9, -1e9], [-np.inf, -2e9], [2e9, np.inf]),  # a step too small to change x
+    ],
+)
+def test_gradient_is_scipys(x, lower, upper):
+    from scipy.optimize._numdiff import approx_derivative
+
+    from cpm.optimisation.fmin import gradient_by_forward_differences
+
+    def f(z):
+        return float(np.sin(z[0]) * z[1] + z[1] ** 2 / 3)
+
+    x, lower, upper = (np.asarray(v, dtype=float) for v in (x, lower, upper))
+    expected = approx_derivative(f, x, method="2-point", abs_step=1e-8, f0=f(x), bounds=(lower, upper))
+    np.testing.assert_array_equal(gradient_by_forward_differences(f, x, f(x), lower, upper, 1e-8), expected)
+
+
+def test_gradient_steps_away_from_a_bound():
+    from cpm.optimisation.fmin import gradient_by_forward_differences
+
+    calls = []
+
+    def f(x):
+        calls.append(x.copy())
+        return float(np.sum(x**2))
+
+    lower, upper = np.array([0.0, 0.0]), np.array([1.0, 1.0])
+    x = np.array([1.0, 0.5])
+    gradient = gradient_by_forward_differences(f, x, f(x), lower, upper, 1e-8)
+    assert all(np.all((c >= lower) & (c <= upper)) for c in calls)
+    assert calls[1][0] < 1.0  # stepped down from the upper bound
+    np.testing.assert_allclose(gradient, 2 * x, rtol=1e-6)
diff --git a/test/optimisation/test_fminbound_hessian.py b/test/optimisation/test_fminbound_hessian.py
new file mode 100644
index 00000000..851d60bd
--- /dev/null
+++ b/test/optimisation/test_fminbound_hessian.py
@@ -0,0 +1,184 @@
+"""
+`FminBound(hessian="finite_differences")`: the Hessian at the optimum by finite
+differences that stay within the bounds, computed once per participant.
+"""
+
+import numdifftools as nd
+import numpy as np
+import pytest
+
+import cpm.optimisation.fmin as fmin_module
+from cpm.core.optimisers import numerical_hessian, objective
+from cpm.optimisation import FminBound, minimise
+from test.optimisation.test_fminbound_gradient import bandit, rlrw
+
+LOSS = minimise.LogLikelihood.bernoulli
+A = np.array([[4.0, 1.5, 0.2], [1.5, 0.7, -0.1], [0.2, -0.1, 2.0]])
+
+
+def counted(f):
+    calls = []
+
+    def g(x):
+        calls.append(np.array(x, copy=True))
+        return f(x)
+
+    return g, calls
+
+
+def quadratic(x):
+    return float(0.5 * x @ A @ x + np.sin(1.0))
+
+
+@pytest.mark.parametrize("d", [2, 3])
+def test_interior_hessian_of_a_quadratic(d):
+    from cpm.core.optimisers import finite_difference_hessian
+
+    f, calls = counted(lambda x: quadratic(np.pad(x, (0, 3 - d))))
+    x = np.array([0.3, 2.0, -1.0])[:d]
+    H = finite_difference_hessian(f, x, np.full(d, -10.0), np.full(d, 10.0))
+    np.testing.assert_allclose(H, A[:d, :d], rtol=1e-4, atol=1e-6)
+    assert len(calls) == d * d + d + 1
+
+
+def test_hessian_on_an_upper_bound_stays_within_the_bounds():
+    from cpm.core.optimisers import finite_difference_hessian
+
+    f, calls = counted(quadratic)
+    lower, upper = np.array([-10.0, 0.0, -10.0]), np.array([10.0, 2.0, 10.0])
+    x = np.array([0.3, 2.0, -1.0])
+    H = finite_difference_hessian(f, x, lower, upper)
+    assert all(np.all((c >= lower) & (c <= upper)) for c in calls)
+    assert len(calls) == 3 * 4 // 2 + 3 + 1
+    np.testing.assert_allclose(H, A, rtol=1e-4, atol=1e-6)
+
+
+def learners(participants=6, trials=120, alpha=0.3, temperature=4.0, seed=2):
+    """A bandit with rewarding arms and the choices of RLRW, so that the estimates are inside the bounds."""
+    rng = np.random.default_rng(seed)
+    data = bandit(participants, trials, seed)
+    reward_probability = np.array([0.2, 0.4, 0.6, 0.8])
+    for side in ("left", "right"):
+        data[f"reward_{side}"] = (rng.random(len(data)) < reward_probability[data[f"arm_{side}"] - 1]).astype(float)
+    choices = []
+    for _, g in data.groupby("ppt"):
+        values = np.full(4, 0.25)
+        for left, right, reward_left, reward_right in g[["arm_left", "arm_right", "reward_left", "reward_right"]].to_numpy():
+            left, right = int(left) - 1, int(right) - 1
+            p_right = 1 / (1 + np.exp(-temperature * (values[right] - values[left])))
+            choice = int(rng.random() < p_right)
+            arm, reward = (right, reward_right) if choice else (left, reward_left)
+            values[arm] += alpha * (reward - values[arm])
+            choices.append(choice)
+    data["response"] = data["observed"] = choices
+    return data
+
+
+def test_hessian_of_rlrw_agrees_with_numdifftools_central():
+    from cpm.core.optimisers import finite_difference_hessian
+
+    data = learners()
+    model = rlrw(data)
+    fit = FminBound(model=model, data=data.groupby("ppt"), minimisation=LOSS, prior=True, number_of_starts=1,
+                    initial_guess=[0.5, 5.0], approx_grad=True)
+    fit.optimise()
+    lower, upper = (np.asarray(b, dtype=float) for b in model.parameters.bounds())
+    checked = 0
+    for (_, participant), result in zip(data.groupby("ppt"), fit.fit):
+        x = result["x"]
+        ## numdifftools' own steps must stay within the bounds too
+        if np.any(x <= lower + 0.05 * (upper - lower)) or np.any(x >= upper - 0.05 * (upper - lower)):
+            continue
+        model.reset(data=participant)
+        observed = participant.observed.to_numpy()
+        f = lambda z: objective(z, model, observed, LOSS, True)  # noqa: E731
+        expected = nd.Hessian(f, method="central")(x)
+        H = finite_difference_hessian(f, x, lower, upper)
+        assert np.abs(H - expected).max() <= 1e-4 * np.abs(expected).max()
+        checked += 1
+    assert checked > 0
+
+
+def fit(hessian, temperature_upper=10.0, starts=2, **kwargs):
+    data = bandit()
+    np.random.seed(3)
+    optimiser = FminBound(
+        model=rlrw(data, temperature_upper),
+        data=data.groupby("ppt"),
+        minimisation=LOSS,
+        prior=True,
+        number_of_starts=starts,
+        approx_grad=True,
+        hessian=hessian,
+        **kwargs,
+    )
+    optimiser.optimise()
+    return optimiser, data
+
+
+def test_default_is_the_numdifftools_hessian_next_to_the_optimum():
+    optimiser, data = fit("numdifftools")
+    model = rlrw(data)
+    for (_, participant), result in zip(data.groupby("ppt"), optimiser.fit):
+        model.reset(data=participant)
+        observed = participant.observed.to_numpy()
+        expected = numerical_hessian(lambda z: objective(z, model, observed, LOSS, True), result["x"] + 1e-3)
+        np.testing.assert_array_equal(result["hessian"], expected)
+
+
+def test_only_the_hessian_changes():
+    old, _ = fit("numdifftools")
+    new, _ = fit("finite_differences")
+    assert len(old.fit) == len(new.fit)
+    for a, b in zip(old.fit, new.fit):
+        assert set(a) == set(b)
+        for key in a:
+            if key != "hessian":
+                np.testing.assert_array_equal(a[key], b[key])
+    assert old.parameters == new.parameters
+
+
+def test_hessian_and_evaluate_fit_once_per_participant(monkeypatch):
+    hessians, evaluations = [], []
+    original_hessian, original_evaluate = fmin_module.finite_difference_hessian, fmin_module.evaluate_fit
+
+    def hessian(*args, **kwargs):
+        hessians.append(1)
+        return original_hessian(*args, **kwargs)
+
+    def evaluate(*args, **kwargs):
+        evaluations.append(1)
+        return original_evaluate(*args, **kwargs)
+
+    monkeypatch.setattr(fmin_module, "finite_difference_hessian", hessian)
+    monkeypatch.setattr(fmin_module, "evaluate_fit", evaluate)
+    optimiser, data = fit("finite_differences", starts=3)
+    participants = data.ppt.nunique()
+    assert len(hessians) == participants
+    assert len(evaluations) == participants
+    assert all(np.all(np.isfinite(r["hessian"])) for r in optimiser.fit)
+
+
+def test_only_the_hessian_changes_in_parallel():
+    old, _ = fit("numdifftools", parallel=True, cl=2)
+    new, _ = fit("finite_differences", parallel=True, cl=2)
+    for a, b in zip(old.fit, new.fit):
+        np.testing.assert_array_equal(a["x"], b["x"])
+        assert a["fun"] == b["fun"]
+        assert np.all(np.isfinite(b["hessian"]))
+
+
+def test_no_all_zero_hessian_on_an_upper_bound():
+    old, _ = fit("numdifftools", temperature_upper=0.5)
+    new, _ = fit("finite_differences", temperature_upper=0.5)
+    on_bound = [i for i, r in enumerate(new.fit) if r["x"][1] == 0.5]
+    assert on_bound
+    assert any(np.all(old.fit[i]["hessian"] == 0) for i in on_bound)  # what this fixes
+    for i in on_bound:
+        assert not np.all(new.fit[i]["hessian"] == 0)
+        np.linalg.cholesky(new.fit[i]["hessian"])
+
+
+def test_unknown_hessian_is_an_error():
+    with pytest.raises(ValueError, match="hessian"):
+        fit("central")
diff --git a/test/optimisation/test_hessian_option.py b/test/optimisation/test_hessian_option.py
new file mode 100644
index 00000000..cecfdec2
--- /dev/null
+++ b/test/optimisation/test_hessian_option.py
@@ -0,0 +1,116 @@
+"""
+The `hessian` option of every optimiser: "numdifftools" (the default, the previous
+Hessian) or "finite_differences" (at the optimum, within the bounds, once per
+participant for the best start).
+"""
+
+import numpy as np
+import pytest
+
+import cpm.optimisation.fmin as fmin_module
+import cpm.optimisation.free as free_module
+from cpm.core.optimisers import finite_difference_hessian, objective
+from cpm.optimisation import Fmin, FminBound, Minimize, minimise
+from test.optimisation.test_fminbound_gradient import bandit, rlrw
+
+try:
+    import cpm.optimisation.bads as bads_module
+    from cpm.optimisation import Bads
+
+    HAS_BADS = True
+except ImportError:  # pragma: no cover - pybads is an optional dependency
+    HAS_BADS = False
+
+LOSS = minimise.LogLikelihood.bernoulli
+OPTIMISERS = [Fmin, FminBound, Minimize]
+MODULES = {Fmin: fmin_module, FminBound: fmin_module, Minimize: free_module}
+if HAS_BADS:
+    OPTIMISERS.append(Bads)
+    MODULES[Bads] = bads_module
+IDS = [cls.__name__ for cls in OPTIMISERS]
+PARTICIPANTS = 3
+## Bads records its own run time
+VARYING = {"total_time"}
+
+
+def build(cls, hessian=None):
+    data = bandit(participants=PARTICIPANTS, trials=30)
+    kwargs = {} if hessian is None else {"hessian": hessian}
+    if cls is FminBound:
+        kwargs["approx_grad"] = True
+    if HAS_BADS and cls is Bads:
+        kwargs["options"] = {"display": "off"}
+    np.random.seed(4)
+    optimiser = cls(model=rlrw(data), data=data.groupby("ppt"), minimisation=LOSS, prior=True, number_of_starts=2, **kwargs)
+    optimiser.optimise()
+    return optimiser, data
+
+
+def records(optimiser):
+    """The per-participant records that hold the Hessian: `fit`, except for `Minimize`, which keeps it in `details`."""
+    return optimiser.details if isinstance(optimiser, Minimize) else optimiser.fit
+
+
+def assert_same(a, b, skip=()):
+    assert set(a) == set(b)
+    for key in set(a) - set(skip) - VARYING:
+        np.testing.assert_array_equal(np.asarray(a[key], dtype=object), np.asarray(b[key], dtype=object), err_msg=key)
+
+
+@pytest.mark.parametrize("cls", OPTIMISERS, ids=IDS)
+def test_default_is_numdifftools(cls):
+    default, _ = build(cls)
+    explicit, _ = build(cls, "numdifftools")
+    assert default.hessian == "numdifftools"
+    for a, b in zip(default.fit, explicit.fit):
+        assert_same(a, b)
+
+
+@pytest.mark.parametrize("cls", OPTIMISERS, ids=IDS)
+def test_only_the_hessian_changes(cls):
+    old, _ = build(cls, "numdifftools")
+    new, _ = build(cls, "finite_differences")
+    for a, b in zip(old.fit, new.fit):
+        assert_same(a, b, skip={"hessian"})
+    for a, b in zip(records(old), records(new)):
+        assert not np.array_equal(a["hessian"], b["hessian"])
+    assert old.parameters == new.parameters
+
+
+@pytest.mark.parametrize("cls", OPTIMISERS, ids=IDS)
+def test_hessian_is_the_finite_difference_one_at_the_optimum(cls):
+    fit, data = build(cls, "finite_differences")
+    model = rlrw(data)
+    lower, upper = (np.asarray(b, dtype=float) for b in model.parameters.bounds())
+    for (_, participant), record in zip(data.groupby("ppt"), records(fit)):
+        model.reset(data=participant)
+        observed = participant.observed.to_numpy()
+        x = record["xopt"] if "xopt" in record else record["x"]  # Fmin records its estimates as xopt
+        expected = finite_difference_hessian(lambda z: objective(z, model, observed, LOSS, True), x, lower, upper)
+        np.testing.assert_array_equal(record["hessian"], expected)
+
+
+@pytest.mark.parametrize("cls", OPTIMISERS, ids=IDS)
+def test_hessian_and_extras_once_per_participant(cls, monkeypatch):
+    module = MODULES[cls]
+    calls = {"hessian": 0, "evaluate_fit": 0}
+    original_hessian, original_evaluate = module.finite_difference_hessian, module.evaluate_fit
+
+    def hessian(*args, **kwargs):
+        calls["hessian"] += 1
+        return original_hessian(*args, **kwargs)
+
+    def evaluate(*args, **kwargs):
+        calls["evaluate_fit"] += 1
+        return original_evaluate(*args, **kwargs)
+
+    monkeypatch.setattr(module, "finite_difference_hessian", hessian)
+    monkeypatch.setattr(module, "evaluate_fit", evaluate)
+    build(cls, "finite_differences")
+    assert calls == {"hessian": PARTICIPANTS, "evaluate_fit": PARTICIPANTS}
+
+
+@pytest.mark.parametrize("cls", OPTIMISERS, ids=IDS)
+def test_unknown_hessian_is_an_error(cls):
+    with pytest.raises(ValueError, match="hessian"):
+        build(cls, "central")
diff --git a/test/optimisation/test_minimise.py b/test/optimisation/test_minimise.py
index d7178b04..f8ea0f8b 100644
--- a/test/optimisation/test_minimise.py
+++ b/test/optimisation/test_minimise.py
@@ -370,6 +370,16 @@ def test_bernoulli_and_continuous_equal_scipy():
     assert LogLikelihood.continuous(predicted=x, observed=y) == -np.sum(norm.logpdf(x, y, 1))
 
 
+def test_bernoulli_equals_scipy_when_every_term_is_finite():
+    from scipy.stats import bernoulli
+
+    rng = np.random.default_rng(6)
+    p = rng.uniform(0, 1, (120, 1))
+    k = rng.integers(0, 2, 120)
+    expected = -np.sum(np.nan_to_num(bernoulli.logpmf(k, np.clip(p.ravel(), 1e-10, 1 - 1e-10)), neginf=-1e100))
+    assert LogLikelihood.bernoulli(predicted=p, observed=k) == expected
+
+
 def test_losses_still_name_nan_and_inf():
     with pytest.raises(ValueError, match="Predicted values contain NaN"):
         LogLikelihood.bernoulli(predicted=np.array([0.5, np.nan]), observed=np.array([1, 0]))