Example 31: excursion set with BSS-style SMC search

Script: examples/example31_excursionset_smc.py

Purpose

The script estimates an excursion set with ExcursionSetBSS. It replaces the fixed grid of example30 by a particle population. The particle target is controlled by an interpolation parameter mu between an initial threshold and the final excursion threshold. The construction is related to Bayesian subset simulation [1].

What is computed

• posterior mean and variance at particle positions.

• the threshold u(mu) = (1 - mu) * u_init + mu * u_target.

• a target log-density proportional to log P(Y(x) > u(mu) | observations) inside the input box.

• SMC reweighting, resampling, and Markov moves of the particle population.

• excursion probabilities and excursion-set criteria at the current threshold.

Main objects

• gpmpcontrib.optim.excursionset.ExcursionSetBSS

• gpmpcontrib.SequentialStrategyBSS

• gpmp.mcmc.smc.SMC

Outputs

Run python examples/example31_excursionset_smc.py from the repository root to execute the example. Regenerate the static figure with cd docs && python make_example_results.py.

Intermediate threshold with mu = 0.85. Top panel: current GP posterior and observations. Middle panel: excursion-set criterion at the current threshold. Lower panel: excursion probability with SMC particles. Particle heights show the current target density up to a multiplicative constant.

Source excerpt

        print("(m[n]) Move particles without updating threshold [n] times (default: 1)")
        print("(u[n]) Move particles and update threshold [n] times (default: 1)")
        print("--")
        print("(e[n]) Make a new evaluation [n] times (default: 1)")
        print("--")
        print("(r) Restart the particles")
        print("(p) Plot current state")
        print("(q) Quit")

        user_input = input("Enter your choice: ").strip().lower()

        # Extract command and optional number or value (supports "m3", "u2", "s 0.4", "n 0.2", etc.)
        match = re.match(r"([muerpsqn])\s*([\d\.]*)", user_input)

        if not match:
            print("Invalid input. Please enter a valid choice.")
            continue

        user_choice, param_str = match.groups()

        # Handle the "s" option separately (setting mu manually)
        if user_choice == "s":
            try:
                mu_value = float(param_str)
                if not (0 <= mu_value <= 1):
                    print("mu must be between 0 and 1.")
                    continue
                algo.mu = mu_value
                print(f"Set mu to {algo.mu:.4f}")
            except ValueError:
                print("Invalid value for mu. Please enter a number between 0 and 1.")
            continue  # Skip to next iteration after setting mu

        # Default repeat count for commands that take an integer (m, u, e)
        repeat = int(param_str) if param_str.isdigit() else 1