# Dynamic occupancy model with pymc\_extras

**URL:** <https://discourse.pymc.io/t/dynamic-occupancy-model-with-pymc-extras/17436>\
**Category:** General\
**Created:** [November 9, 2025, 8:04pm UTC](https://discourse.pymc.io/t/dynamic-occupancy-model-with-pymc-extras/17436 "2025-11-09T20:04:56Z")\
**Posts on this page:** 17\
**Page:** 1

<div class="post-metadata">

**Author:** ![philpatton](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/philpatton/32/7674_2.png) [@philpatton](https://discourse.pymc.io/u/philpatton)\
**Post date:** [November 9, 2025, 8:04pm UTC](https://discourse.pymc.io/t/dynamic-occupancy-model-with-pymc-extras/17436/1 "2025-11-09T20:04:56Z")

</div>

Good afternoon!

I’m writing to see if y’all think that this dynamic occupancy model would be a good fit for `DiscreteMarkovChain` and `marginalize` in the PyMC Extras package. In this model, we detect animals at occupied, z=1, sites (aka patches) with probability p. We assume that animals can’t be detected at unoccupied, z=0, sites. Between seasons, unoccupied sites can be colonized with probability \gamma, and occupied sites can go extinct with probability \epsilon. Within a season, sites are surveyed and closed to colonization / extinction dynamics. Below is some simulation code, and an example of how to fit the model in NumPyro. While the NumPyro version works great, I’m curious if this would be a good fit for the goodies in pymc-extras. Personally I find PyMC syntax a little more digestible.

Thanks in advance for your help!  
Phil

```auto
from jax import random
from numpyro.contrib.control_flow import scan
from numpyro.infer import NUTS, MCMC, Predictive
import arviz as az
import jax.numpy as jnp
import numpy as np
import numpyro
import numpyro.distributions as dist

# hyperparameters
RANDOM_SEED = 1792

## true values for colext model
PSI_TRUE = 0.6
EPSILON_TRUE = 0.3
GAMMA_TRUE = 0.15
P_TRUE = 0.4
SITE_COUNT = 250
SURVEY_COUNT = 3
SEASON_COUNT = 10
interval_count = SEASON_COUNT - 1

def simulate_data():
    """Simulate detection/non-detection data from a dynamic occupancy model"""

    rng = np.random.default_rng(RANDOM_SEED)

    # empty array to fill in the occupancy states later
    z = np.zeros((SITE_COUNT, SEASON_COUNT), dtype=int)

    # initial values for the occupancy state
    z[:, 0] = rng.binomial(n=1, p=PSI_TRUE, size=SITE_COUNT)

    # simulate transitions
    for t in range(1, SEASON_COUNT):

        # patches can be colonized, go extinct, remain occupied, or remain unoccupied
        mu_z = z[:, t-1] * (1 - EPSILON_TRUE) + (1 - z[:, t-1]) * GAMMA_TRUE
        z[:, t] = rng.binomial(n=1, p=mu_z)

    # simulate detection non-detection data
    mu_x = z * P_TRUE
    x = rng.binomial(n=1, p=mu_x[:, :, None],
                    size=(SITE_COUNT, SEASON_COUNT, SURVEY_COUNT))
    return x

def dynamic_occupancy(detection_history):
    '''Dynamic occupancy model in NumPyro.'''
    site_count, season_count, survey_count = detection_history.shape

    # scalar priors for the four probabilistic parameters
    psi = numpyro.sample("psi", dist.Uniform(0, 1)) # initial occupancy prob
    gamma = numpyro.sample("gamma", dist.Uniform(0, 1)) # colonization prob
    epsilon = numpyro.sample("epsilon", dist.Uniform(0, 1)) # extinction prob
    p = numpyro.sample("p", dist.Uniform(0, 1)) # recapture prob

    def transition_and_detect(carry, y_t):
        """Transitions betweens states and defines the likelihood."""

        # unpack the values that are returned from the transition function at
        # the previous time step
        z_prev, t = carry

        # transition the latent state at every site
        with numpyro.plate("sites", site_count):

            # probability of transitioning according to the previous state
            mu_z_t = z_prev * (1 - epsilon) + (1 - z_prev) * gamma

            # dist.util.clamp_probs() helps the sampler avoid boundary regions
            z = numpyro.sample(
                "z",
                dist.Bernoulli(dist.util.clamp_probs(mu_z_t)),
                infer={"enumerate": "parallel"}, # this is where we marginalize!
            )

            # the likelihood of each observation at each site
            mu_y = z * p
            with numpyro.plate('surveys', survey_count):
                numpyro.sample(
                    "y",
                    dist.Bernoulli(dist.util.clamp_probs(mu_y)),
                    obs=y_t.T
                )

            # carry forward the current z state and incremented time index
            # None indicates we don't return/accumulate any outputs from scan
            return (z, t + 1), None

    # the initial state only depends on psi
    with numpyro.plate('sites', site_count):
        z0 = numpyro.sample(
            "z0",
            dist.Bernoulli(dist.util.clamp_probs(psi)),
            infer={"enumerate": "parallel"},
        )

        # compute the likelihood of the detection data for just the first season
        mu_y = z0 * p
        with numpyro.plate('surveys', survey_count):
            numpyro.sample(
                "y0",
                dist.Bernoulli(dist.util.clamp_probs(mu_y)),
                obs=detection_history[:, 0].T # just the first occasion!
            )

    # now we scan (or apply) the transition function across the remaining seasons
    scan(
        transition_and_detect, # function to scan
        (z0, 0), # initial states
        jnp.swapaxes(detection_history[:, 1:], 0, 1), # scan across first dimension of data
    )

detections = simulate_data()
rng_key = random.PRNGKey(RANDOM_SEED)

# specify which sampler you want to use
nuts_kernel = NUTS(dynamic_occupancy)

# configure the MCMC run
mcmc = MCMC(nuts_kernel, num_warmup=500, num_samples=1000, num_chains=4)

# run the MCMC then inspect the output
mcmc.run(rng_key, detections)
mcmc.print_summary()

```

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [November 9, 2025, 8:39pm UTC](https://discourse.pymc.io/t/dynamic-occupancy-model-with-pymc-extras/17436/2 "2025-11-09T20:39:09Z")

</div>

I think it would work

---

<div class="post-metadata">

**Author:** ![jessegrabowski](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/jessegrabowski/32/5010_2.png) [@jessegrabowski](https://discourse.pymc.io/u/jessegrabowski)\
**Post date:** [November 9, 2025, 11:02pm UTC](https://discourse.pymc.io/t/dynamic-occupancy-model-with-pymc-extras/17436/3 "2025-11-09T23:02:01Z")

</div>

It would make a good case study for pymc-examples 🙂

---

<div class="post-metadata">

**Author:** ![philpatton](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/philpatton/32/7674_2.png) [@philpatton](https://discourse.pymc.io/u/philpatton)\
**Post date:** [November 9, 2025, 11:48pm UTC](https://discourse.pymc.io/t/dynamic-occupancy-model-with-pymc-extras/17436/4 "2025-11-09T23:48:35Z")

</div>

Awesome! Well, funnily enough, the act of asking the question seemed to solve my problem! I wasn’t able to get a previous version the code below to work. But it turns out I was just missing the shape argument. Now it’s working and returns parameters.

I would be happy to add this case study to pymc-examples if that would be helpful (I actually have several other ecological examples in PyMC handy). Sadly I know nothing about the process. What does it entail? Do I just submit a Jupyter Notebook as a pull request to pymc-examples?

```auto
with pm.Model() as colext:

    # simple model with constant probabilities
    psi = pm.Uniform('ψ', 0, 1)
    epsilon = pm.Uniform('ε', 0, 1)
    gamma = pm.Uniform('γ', 0, 1)
    p = pm.Uniform('p', 0, 1)

    transition_matrix = pt.stack(
            [[1 - gamma, gamma],
             [epsilon, 1 - epsilon]], axis=0
    )

    # initial occupancy only depends on psi
    initial_distribution = pm.Bernoulli.dist(psi, shape=SITE_COUNT)

    Z = pmx.DiscreteMarkovChain(
        'Z',
        transition_matrix,
        steps=interval_count,
        init_dist=initial_distribution,
        shape=(SITE_COUNT, SEASON_COUNT)
    )

    mu_x = Z[:, :, None] * p
    pm.Bernoulli('x', mu_x, observed=x)

colext_marginal = pmx.marginalize(colext, ['Z'])

with colext_marginal:
    idata = pm.sample()

```

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [November 10, 2025, 8:45am UTC](https://discourse.pymc.io/t/dynamic-occupancy-model-with-pymc-extras/17436/5 "2025-11-10T08:45:53Z")

</div>

That’s quite readable in the end.

To publish to pymc-examples it is suggested you first open a PR that gives some reviewers context for the Notebook. There’s a specific template for new notebook propoosals: [GitHub · Where software is built](https://github.com/pymc-devs/pymc-examples/issues)

Then open a PR following these guidelines (somethings are unrelated to you, just skim): [GitHub - pymc-devs/pymc-examples: Examples of PyMC models, including a library of Jupyter notebooks.](https://github.com/pymc-devs/pymc-examples/?tab=contributing-ov-file)

You can tag me or @jessegrabowski for review. There are some formatting guidelines / gotchas but nothing too cumbersome (I hope).

Ecology examples would be superb. Just check there’s nothing too similar already. If there is we can update/merge/replace if your material is better.

---

<div class="post-metadata">

**Author:** ![philpatton](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/philpatton/32/7674_2.png) [@philpatton](https://discourse.pymc.io/u/philpatton)\
**Post date:** [July 9, 2026, 5:18pm UTC](https://discourse.pymc.io/t/dynamic-occupancy-model-with-pymc-extras/17436/6 "2026-07-09T17:18:36Z")

</div>

Hey Y’all! Super excited to see the time-varying parameters added to DiscreteMarkovChain! This opens up a lot of possibilities for classical ecological models (e.g., capture-recapture).

I know that the marginalize portion of PyMC-extras is undergoing some change, but I just wanted to flag anyway that the code above no longer seems to work. I tried a fresh re-install of pymc (version `6.1.0`) and pymc-extras (version `0.12.0`) and got this error message:

```auto
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
Cell In[1], line 76
     72 
     73 colext_marginal = pmx.marginalize(colext, ['Z'])
     74 
     75 with colext_marginal:
---> 76 idata = pm.sample()
     77 
     78 az.plot_trace(idata)

File ~/source/repos/popan/.pixi/envs/default/lib/python3.14/site-packages/pymc/sampling/mcmc.py:917, in sample(draws, tune, chains, cores, random_seed, progressbar, progressbar_theme, quiet, step, var_names, nuts_sampler, initvals, init, jitter_max_retries, n_init, trace, discard_tuned_samples, compute_convergence_checks, keep_warning_stat, return_inferencedata, idata_kwargs, nuts_sampler_kwargs, callback, mp_ctx, blas_cores, model, backend, compile_kwargs, **kwargs)
    914 msg = f"Only {draws} samples per chain. Reliable r-hat and ESS diagnostics require longer chains for accurate estimate."
    915 _log.warning(msg)
--> 917 provided_steps, selected_steps = assign_step_methods(model, step, methods=STEP_METHODS)
    918 exclusive_nuts = (
    919 # User provided an instantiated NUTS step, and nothing else is needed
    920 (not selected_steps and len(provided_steps) == 1 and isinstance(provided_steps[0], NUTS))
   (...) 927 )
    928 )
    930 if nuts_sampler is None:
    931 # Try to use nutpie by default if no setting is clearly at odds.
    932 # Requires all model variables, numba or jax preference,
    933 # and must not conflict with pymc sample-only arguments.

File ~/source/repos/popan/.pixi/envs/default/lib/python3.14/site-packages/pymc/sampling/mcmc.py:324, in assign_step_methods(model, step, methods)
    322 methods_list: list[type[BlockedStep]] = list(methods or STEP_METHODS)
    323 selected_steps: dict[type[BlockedStep], list] = {}
--> 324 model_logp = model.logp()
    326 for var in model.value_vars:
    327 if var not in assigned_vars:
    328 # determine if a gradient can be computed

File ~/source/repos/popan/.pixi/envs/default/lib/python3.14/site-packages/pymc/model/core.py:728, in Model.logp(self, vars, jacobian, sum)
    726 rv_logps: list[TensorVariable] = []
    727 if rvs:
--> 728 rv_logps = transformed_conditional_logp(
    729 rvs=rvs,
    730 rvs_to_values=self.rvs_to_values,
    731 rvs_to_transforms=self.rvs_to_transforms,
    732 jacobian=jacobian,
    733 )
    734 assert isinstance(rv_logps, list)
    736 # Replace random variables by their value variables in potential terms

File ~/source/repos/popan/.pixi/envs/default/lib/python3.14/site-packages/pymc/logprob/basic.py:642, in transformed_conditional_logp(rvs, rvs_to_values, rvs_to_transforms, jacobian, **kwargs)
    639 transform_rewrite = TransformValuesRewrite(values_to_transforms) # type: ignore[arg-type]
    641 kwargs.setdefault("warn_rvs", False)
--> 642 temp_logp_terms = conditional_logp(
    643 rvs_to_values,
    644 extra_rewrites=transform_rewrite,
    645 use_jacobian=jacobian,
    646 **kwargs,
    647 )
    649 # The function returns the logp for every single value term we provided to it.
    650 # This includes the extra values we plugged in above, so we filter those we
    651 # actually wanted in the same order they were given in.
    652 logp_terms = {}

File ~/source/repos/popan/.pixi/envs/default/lib/python3.14/site-packages/pymc/logprob/basic.py:572, in conditional_logp(rv_values, warn_rvs, ir_rewriter, extra_rewrites, **kwargs)
    569 node_values = remapped_vars[: len(node_values)]
    570 node_inputs = remapped_vars[len(node_values) :]
--> 572 node_logprobs = _logprob(
    573 node.op,
    574 node_values,
    575 *node_inputs,
    576 **kwargs,
    577 )
    579 if not isinstance(node_logprobs, list | tuple):
    580 node_logprobs = [node_logprobs]

File ~/source/repos/popan/.pixi/envs/default/lib/python3.14/functools.py:982, in singledispatch.<locals>.wrapper(*args, **kw)
    979 if not args:
    980 raise TypeError(f'{funcname} requires at least '
    981 '1 positional argument')
--> 982 return dispatch(args[0]. __class__ )(*args, **kw)

File ~/source/repos/popan/.pixi/envs/default/lib/python3.14/site-packages/pymc_extras/model/marginal/distributions.py:308, in marginal_hmm_logp(op, values, *inputs, **kwargs)
    306 @_logprob.register(MarginalDiscreteMarkovChainRV)
    307 def marginal_hmm_logp(op, values, *inputs, **kwargs):
--> 308 chain_rv, *dependent_rvs = inline_ofg_outputs(op, inputs)
    310 P, n_steps_, init_dist_, rng = chain_rv.owner.inputs
    311 domain = pt.arange(P.shape[-1], dtype="int32")

File ~/source/repos/popan/.pixi/envs/default/lib/python3.14/site-packages/pymc_extras/model/marginal/distributions.py:217, in inline_ofg_outputs(op, inputs)
    211 def inline_ofg_outputs(op: OpFromGraph, inputs: Sequence[Variable]) -> tuple[Variable]:
    212 """Inline the inner graph (outputs) of an OpFromGraph Op.
    213 
    214 Whereas `OpFromGraph` "wraps" a graph inside a single Op, this function "unwraps"
    215 the inner graph.
    216 """
--> 217 return graph_replace(
    218 op.inner_outputs,
    219 replace=tuple(zip(op.inner_inputs, inputs)),
    220 strict=False,
    221 )

File ~/source/repos/popan/.pixi/envs/default/lib/python3.14/site-packages/pytensor/graph/replace.py:160, in graph_replace(outputs, replace, strict)
    156 equiv = {c: c.clone(name=f"i-{i}") for i, c in enumerate(conditions)}
    157 # some replace keys may disappear
    158 # the reason is they are outside the graph
    159 # clone the graph but preserve the equiv mapping
--> 160 fg = FunctionGraph(
    161 conditions,
    162 outputs,
    163 # clone_get_equiv kwargs
    164 copy_orphans=False,
    165 copy_inputs=False,
    166 memo=equiv,
    167 )
    168 # replace the conditions back
    169 fg_replace = {equiv[c]: c for c in conditions}

File ~/source/repos/popan/.pixi/envs/default/lib/python3.14/site-packages/pytensor/graph/fg.py:139, in FunctionGraph. __init__ (self, inputs, outputs, features, clone, update_mapping, **clone_kwds)
    134 inputs = [
    135 i for i in graph_inputs(outputs) if not isinstance(i, AtomicVariable)
    136 ]
    138 if clone:
--> 139 _memo = clone_get_equiv(
    140 inputs,
    141 outputs,
    142 **clone_kwds,
    143 )
    144 outputs = [cast(Variable, _memo[o]) for o in outputs]
    145 inputs = [cast(Variable, _memo[i]) for i in inputs]

File ~/source/repos/popan/.pixi/envs/default/lib/python3.14/site-packages/pytensor/graph/basic.py:1057, in clone_get_equiv(inputs, outputs, copy_inputs, copy_orphans, memo, clone_inner_graphs, **kwargs)
   1054 else:
   1055 memo[input] = input
-> 1057 clone_node_and_cache(apply, memo, **kwargs)
   1059 # finish up by cloning any remaining outputs (it can happen)
   1060 for output in outputs:

File ~/source/repos/popan/.pixi/envs/default/lib/python3.14/site-packages/pytensor/graph/basic.py:980, in clone_node_and_cache(node, clone_d, **kwargs)
    976 return None
    978 cloned_inputs: list[Variable] = [cast(Variable, clone_d[i]) for i in node.inputs]
--> 980 new_node = node.clone_with_new_inputs(cloned_inputs, **kwargs)
    982 clone_d[node] = new_node
    984 for old_o, new_o in zip(node.outputs, new_node.outputs, strict=True):

AttributeError: 'FrozenApply' object has no attribute 'clone_with_new_inputs'

```

Here was the code block.

```auto
import arviz as az
import numpy as np
import pymc as pm
import pymc_extras as pmx 
import pytensor.tensor as pt

# hyperparameters
RANDOM_SEED = 1792

## true values for colext model
PSI_TRUE = 0.6
EPSILON_TRUE = 0.3
GAMMA_TRUE = 0.15
P_TRUE = 0.4
SITE_COUNT = 250
SURVEY_COUNT = 3
SEASON_COUNT = 10
interval_count = SEASON_COUNT - 1

def simulate_data():
    """Simulate detection/non-detection data from a dynamic occupancy model"""

    rng = np.random.default_rng(RANDOM_SEED)

    # empty array to fill in the occupancy states later
    z = np.zeros((SITE_COUNT, SEASON_COUNT), dtype=int)

    # initial values for the occupancy state
    z[:, 0] = rng.binomial(n=1, p=PSI_TRUE, size=SITE_COUNT)

    # simulate transitions
    for t in range(1, SEASON_COUNT):

        # patches can be colonized, go extinct, remain occupied, or remain unoccupied
        mu_z = z[:, t-1] * (1 - EPSILON_TRUE) + (1 - z[:, t-1]) * GAMMA_TRUE
        z[:, t] = rng.binomial(n=1, p=mu_z)

    # simulate detection non-detection data
    mu_x = z * P_TRUE
    x = rng.binomial(n=1, p=mu_x[:, :, None],
                    size=(SITE_COUNT, SEASON_COUNT, SURVEY_COUNT))
    return x

x = simulate_data()

with pm.Model() as colext:

    # simple model with constant probabilities
    psi = pm.Uniform('ψ', 0, 1)
    epsilon = pm.Uniform('ε', 0, 1)
    gamma = pm.Uniform('γ', 0, 1)
    p = pm.Uniform('p', 0, 1)

    transition_matrix = pt.stack(
            [[1 - gamma, gamma],
             [epsilon, 1 - epsilon]], axis=0
    )

    # initial occupancy only depends on psi
    initial_distribution = pm.Bernoulli.dist(psi, shape=SITE_COUNT)

    Z = pmx.DiscreteMarkovChain(
        'Z',
        transition_matrix,
        steps=interval_count,
        init_dist=initial_distribution,
        shape=(SITE_COUNT, SEASON_COUNT)
    )

    mu_x = Z[:, :, None] * p
    pm.Bernoulli('x', mu_x, observed=x)

colext_marginal = pmx.marginalize(colext, ['Z'])

with colext_marginal:
    idata = pm.sample()

az.plot_trace(idata)

```

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [July 9, 2026, 6:00pm UTC](https://discourse.pymc.io/t/dynamic-occupancy-model-with-pymc-extras/17436/7 "2026-07-09T18:00:42Z")

</div>

We are in the process of making the dependencies compatible (and pinning it strictly, pymc-extras was too losely pinned, now that’s a bit less experimental we are going to pin better): [Compatibility with PyTensor v3.1 and PyMC v6.1 by ricardoV94 · Pull Request #708 · pymc-devs/pymc-extras · GitHub](https://github.com/pymc-devs/pymc-extras/pull/708)

Also trying to support higher lags: [Support marginalization of HMM with lags\>1 by ricardoV94 · Pull Request #707 · pymc-devs/pymc-extras · GitHub](https://github.com/pymc-devs/pymc-extras/pull/707)

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [July 9, 2026, 6:02pm UTC](https://discourse.pymc.io/t/dynamic-occupancy-model-with-pymc-extras/17436/8 "2026-07-09T18:02:44Z")

</div>

Feedback on the time-varing would be nice. For now we focused only on behavior and haven’t benchmarked / optimized for speed.

I’m curious how people actually want to use time-varying, is it okay to expect a dense vector of transition probabilities, or would a function that creates the transition per iteration (based on the counter / state) be more natural?

The reason we added it was to support conditional distribution after marginalization, which require time-varying Ps

---

<div class="post-metadata">

**Author:** ![philpatton](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/philpatton/32/7674_2.png) [@philpatton](https://discourse.pymc.io/u/philpatton)\
**Post date:** [July 9, 2026, 10:34pm UTC](https://discourse.pymc.io/t/dynamic-occupancy-model-with-pymc-extras/17436/9 "2026-07-09T22:34:49Z")

</div>

Awesome!! I’m glad to hear that it’s all coming together.

That’s a great question. Ecologists like to model the transition probabilities with predictors (e.g., weather), or as time-varying fixed (or random) effects like

```auto
def dynamic_occupancy_model(data):

    ...

    with numpyro.plate('seasons', season_count - 1):
        gamma = numpyro.sample('gamma', dist.Uniform(0, 1))
        epsilon = numpyro.sample('epsilon', dist.Normal(0, 1)) 

    def transition_and_observe(carry, x_current):

        z_previous, t = carry

        # transition probability matrix
        trans_probs = jnp.array([
            [1 - gamma[t], gamma[t]], 
            [epsilon[t], 1 - epsilon[t]]  
        ])

        with numpyro.plate('sites', site_count):

            ...

```

Because there’s spatial replication, they also frequently model the effect of predictors (e.g., forest cover) on transition probabilities. So there might just be a handful of parameters but could end up with `gamma` and `epsilon` having shape `(season_count - 1, site_count)`. Personally I don’t mind stacking them into a big array `P` but perhaps a function would be more intuitive. But I have absolutely zero-skill when it comes to software design so I trust y’all’s intuition!

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [July 10, 2026, 8:11am UTC](https://discourse.pymc.io/t/dynamic-occupancy-model-with-pymc-extras/17436/10 "2026-07-10T08:11:30Z")

</div>

The ultimate design goal is not to force too much into DiscreteMarkovChain but to allow writing arbitrary Scans with categorical emissions and auto-marginalize/recover them. Then users are free to define the transition probabilities outside/inside and pass whatever predictors they want as sequences or loop invariants.

For now that’s just a wish though …

---

<div class="post-metadata">

**Author:** ![philpatton](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/philpatton/32/7674_2.png) [@philpatton](https://discourse.pymc.io/u/philpatton)\
**Post date:** [July 10, 2026, 4:37pm UTC](https://discourse.pymc.io/t/dynamic-occupancy-model-with-pymc-extras/17436/11 "2026-07-10T16:37:08Z")

</div>

Sweet! That sounds like a great goal. That’s basically the way I interact with NumPyro, and why [I’ve found](https://philpatton.github.io/numpyro-js.html) it to be so useful for HMM’s in ecology. But I would be thrilled to get back to PyMC and nutpie, if only to be able to run things on my Mac again!

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [July 13, 2026, 1:48pm UTC](https://discourse.pymc.io/t/dynamic-occupancy-model-with-pymc-extras/17436/12 "2026-07-13T13:48:09Z")

</div>

@philpatton if you want to give it a try, it should now work with the latest releases: [Release v0.13.0 · pymc-devs/pymc-extras · GitHub](https://github.com/pymc-devs/pymc-extras/releases/tag/v0.13.0)

---

<div class="post-metadata">

**Author:** ![philpatton](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/philpatton/32/7674_2.png) [@philpatton](https://discourse.pymc.io/u/philpatton)\
**Post date:** [July 13, 2026, 4:46pm UTC](https://discourse.pymc.io/t/dynamic-occupancy-model-with-pymc-extras/17436/13 "2026-07-13T16:46:32Z")

</div>

Awesome, it works! And I’m thrilled to see [Allow closed form marginalization and implement `conditional` model transform by ricardoV94 · Pull Request #688 · pymc-devs/pymc-extras · GitHub](https://github.com/pymc-devs/pymc-extras/pull/688). This is actually very helpful for capture-recapture, because we often want to recover the discrete alive/dead state [conditional on the capture history](https://forum.pyro.ai/t/recover-discrete-latent-states-after-enumerate-scan/9018). (Basically, if you saw an individual in 2012, 2014, and 2015, you know they must have been alive in 2013, but you are unsure if they were alive in, say, 2011 or 2016.) Thanks again for the help!

---

<div class="post-metadata">

**Author:** ![jessegrabowski](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/jessegrabowski/32/5010_2.png) [@jessegrabowski](https://discourse.pymc.io/u/jessegrabowski)\
**Post date:** [July 13, 2026, 5:27pm UTC](https://discourse.pymc.io/t/dynamic-occupancy-model-with-pymc-extras/17436/14 "2026-07-13T17:27:35Z")

</div>

If you want to do a capture-recapture example for pymc-examples using this feature it’d be very welcome 🙂

---

<div class="post-metadata">

**Author:** ![philpatton](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/philpatton/32/7674_2.png) [@philpatton](https://discourse.pymc.io/u/philpatton)\
**Post date:** [July 14, 2026, 2:01pm UTC](https://discourse.pymc.io/t/dynamic-occupancy-model-with-pymc-extras/17436/15 "2026-07-14T14:01:55Z")

</div>

I’d be thrilled to! I can whip together a [Jolly-Seber example](https://philpatton.github.io/pymc-jssa.html), which will be good for showing off `conditional` because we want the conditional distribution of the total population size over time given the capture histories.

I’d also like to make use of `time_varying_P` since, for better or for worse, this is [canon in capture-recapture](https://num.pyro.ai/en/0.5.0/examples/capture_recapture.html). But I’m having trouble getting the shapes rights with the replicated time series (see below for an example with the occupancy model, which is a nearly identical model). Any ideas?

```auto
import arviz as az
import numpy as np
import pymc as pm
import pytensor.tensor as pt

from pymc_extras.marginal import marginalize, conditional
from pymc_extras import DiscreteMarkovChain

# hyperparameters
RANDOM_SEED = 1792

## true values for colext model
PSI_TRUE = 0.6
EPSILON_TRUE = 0.3
GAMMA_TRUE = 0.15
P_TRUE = 0.4
SITE_COUNT = 250
SURVEY_COUNT = 3
SEASON_COUNT = 10
interval_count = SEASON_COUNT - 1

def simulate_data():
    """Simulate detection/non-detection data from a dynamic occupancy model"""

    rng = np.random.default_rng(RANDOM_SEED)

    # empty array to fill in the occupancy states later
    z = np.zeros((SITE_COUNT, SEASON_COUNT), dtype=int)

    # initial values for the occupancy state
    z[:, 0] = rng.binomial(n=1, p=PSI_TRUE, size=SITE_COUNT)

    # simulate transitions
    for t in range(1, SEASON_COUNT):

        # patches can be colonized, go extinct, remain occupied, or remain unoccupied
        mu_z = z[:, t-1] * (1 - EPSILON_TRUE) + (1 - z[:, t-1]) * GAMMA_TRUE
        z[:, t] = rng.binomial(n=1, p=mu_z)

    # simulate detection non-detection data
    mu_x = z * P_TRUE
    x = rng.binomial(n=SURVEY_COUNT, p=mu_x, size=(SITE_COUNT, SEASON_COUNT))
    return x, z

x, z_true = simulate_data()

with pm.Model() as colext:

    # simple model with constant probabilities
    psi = pm.Uniform('ψ', 0, 1)
    epsilon = pm.Uniform('ε', 0, 1, shape=interval_count)
    gamma = pm.Uniform('γ', 0, 1, shape=interval_count)
    p = pm.Uniform('p', 0, 1)

    # initial occupancy only depends on psi
    init_dist = pm.Bernoulli.dist(psi, shape=SITE_COUNT)

    # create the array of transition matrices by first stacking each row
    P = pt.stack([
        pt.stack([1 - gamma, gamma], axis=-1),
        pt.stack([epsilon, 1 - epsilon], axis=-1)
    ], axis=-2) # resulting shape (interval_count, 2, 2)

    Z = DiscreteMarkovChain(
        'Z',
        P=P,
        steps=interval_count,
        init_dist=init_dist,
        time_varying_P=True,
        shape=(SITE_COUNT, SEASON_COUNT)
    )

    pm.Binomial('x', SURVEY_COUNT, Z * p, observed=x)

colext_marginal = marginalize(colext, ['Z'])

with colext_marginal:
    idata = pm.sample()

```

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [July 14, 2026, 2:57pm UTC](https://discourse.pymc.io/t/dynamic-occupancy-model-with-pymc-extras/17436/16 "2026-07-14T14:57:49Z")

</div>

Ah that’s a shame. I think it’s fixed in: [Support marginalization of HMM with lags\>1 by ricardoV94 · Pull Request #707 · pymc-devs/pymc-extras · GitHub](https://github.com/pymc-devs/pymc-extras/pull/707)

Let me try to clean that up over the finish line

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [July 19, 2026, 8:53pm UTC](https://discourse.pymc.io/t/dynamic-occupancy-model-with-pymc-extras/17436/17 "2026-07-19T20:53:51Z")

</div>

Hopefully this fixed it @philpatton: [Release v0.13.1 · pymc-devs/pymc-extras · GitHub](https://github.com/pymc-devs/pymc-extras/releases/tag/v0.13.1)