# Fitting a value-based decision model with trial-level context variables

**URL:** <https://discuss.bayesflow.org/t/fitting-a-value-based-decision-model-with-trial-level-context-variables/194>\
**Category:** General\
**Created:** [August 27, 2025, 10:06am UTC](https://discuss.bayesflow.org/t/fitting-a-value-based-decision-model-with-trial-level-context-variables/194 "2025-08-27T10:06:35Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![ErikStuchly](https://avatars.discourse-cdn.com/v4/letter/e/59ef9b/32.png) [@ErikStuchly](https://discuss.bayesflow.org/u/ErikStuchly)\
**Post date:** [August 27, 2025, 10:06am UTC](https://discuss.bayesflow.org/t/fitting-a-value-based-decision-model-with-trial-level-context-variables/194/1 "2025-08-27T10:06:35Z")

</div>

Hi there,

I am very new to amortised bayesian inference, so apologies in advance if my questions seem uninformed or contain some misconceptions!

Anyway, I am hoping to create and fit a decision-making model, not quite like the DDM but one that also generates response times and choice variables, and which consists of 6 parameters. So far so good - my task is a two-alternative social decision task though, where the decision difficulty and thus the choice/RT varies based on the current option pair (i.e., in addition to the 6 parameters, it is also determined by the value of the _“selfish” option/self\_option_, the value of the _“altruistic” option/alt\_option_, and a participant-level variable called the _“preferred allocation”/PA_). Based on my understanding, it would make sense to train the neural network over different trials (i.e., different combinations of self\_option, alt\_option and PA), to make sure the network learns “proper” relationships between the parameters and the choices/RTs conditioned over the context. To make matters slightly more complex, the exact combination of self\_opt, alt\_opt and PA differ across participants (meaning that participants do not share the same decision option set/context), so I really need to train the network on all sorts of combinations of these variables. My question is, what is the proper way to do this? It seems to me that these variables (self\_opt, alt\_opt, PA) should be simulated in a separate function and then fed to the simulator object as a context variable. To that end, I have created the following function which creates an (experimentally valid) combination of self\_opt, alt\_opt and PA:

```python
def make_optionpairs(batch_size=1):
    step = 0.05
    max_int = 19

    # For batch sampling
    alt_ints = np.random.randint(0, max_int, batch_size)
    self_ints = alt_ints + 1 + np.random.randint(0, max_int - alt_ints)
    alt_option = alt_ints * step
    self_option = self_ints * step
    PA = np.random.random(batch_size)

    # Best: always return arrays of shape (batch_size,)
    return {
        "alt_option": alt_option, # shape (batch_size,)
        "self_option": self_option,
        "PA": PA
    }

```

Then I add my prior-generating function and the trial-generating function (for the time being, I set it to generate and put out one choice/RT combination at a time):

```python
def prior_sampler():
    # You can change these to realistic ranges
    n_sims = int(np.round(np.random.uniform(1,11)))
    drift_scalings = np.random.uniform(0.1, 10)
    thresholds = np.random.uniform(0.7, 3.0)
    ndts = np.random.uniform(0,0.3)
    rel_sps = np.random.uniform(0.3,0.7)
    rt_multipliers = np.random.uniform(0.1,100)
    return {"n_sims": n_sims, "drift_scaling": drift_scalings, "threshold": thresholds, "ndt": ndts, "rel_sp": rel_sps, "rt_multiplier": rt_multipliers}

def random_dm(alt_option, self_option, PA,n_sims, drift_scaling, threshold, ndt, rel_sp, rt_multiplier, noise_constant=1, dt=0.001, max_rt=10):
    v = (drift_scaling *
         ((1 - (alt_option - PA)**2) -
          (1 - (self_option - PA)**2))
        )
    drift_rates = v
    drifts = np.broadcast_to(drift_rates, (n_sims, 1)).T
    n_trials = drifts.shape[0]
    acc = np.full(drifts.shape, np.nan)
    rt = np.full(drifts.shape, np.nan)
    max_tsteps = int(max_rt / dt)
    
    x = np.ones(drifts.shape) * rel_sp * threshold # evidence at t=0
    ongoing = np.full(drifts.shape, True)
    
    tstep = 0
    while np.any(ongoing) and tstep < max_tsteps:
        # apply drift and noise only to ongoing trials
        moving_indices = np.where(ongoing)
        x[moving_indices] += (
            drifts[moving_indices] * dt +
            np.random.normal(loc=0, scale=noise_constant * np.sqrt(dt), size=len(moving_indices[0]))
        )
        tstep += 1
        
        # handle threshold crossing
        ended_correct = (x >= threshold)
        ended_incorrect = (x <= 0)
        
        # For correct responses
        mask = (ended_correct & ongoing)
        if np.any(mask):
            acc[mask] = 1
            rt[mask] = dt * tstep
            ongoing[mask] = False

        # For incorrect responses
        mask = (ended_incorrect & ongoing)
        if np.any(mask):
            acc[mask] = -1
            rt[mask] = dt * tstep
            ongoing[mask] = False

    rt = rt / rt_multiplier

    # Aggregate "majority" choice per trial (as in rowSums)
    row_acc = np.nansum(acc, axis=1)
    acc_averaged = np.zeros(n_trials, dtype=int)
    acc_averaged[row_acc > 0] = 1
    acc_averaged[row_acc < 0] = 0
    undecided = (row_acc == 0)
    acc_averaged[undecided] = np.random.choice([1, 0], size=np.sum(undecided))

    rt_averaged = np.nansum(rt, axis=1) + ndt

    return dict(accuracy=acc_averaged, rt=rt_averaged)

```

And finally combine it into one simulator

`simulator = bf.make_simulator([prior_sampler, random_dm], meta_fn=make_optionpairs)`

Here are my questions:

1. What I am unsure about is that, when I feed this all into the adapter like so:

```python
adapter = (
    bf.Adapter()
    .to_array()
    #.constrain(["n_sims","drift_scaling", "threshold", "ndt", "rel_sp", "rt_multiplier"],lower=np.array([0,0,0,0,0,0]),upper=np.array([100,10,10,1.5,1,5]))
    .as_set(["accuracy", "rt"])
    .convert_dtype("float64", "float32")
    .concatenate(["accuracy","rt"], into="summary_variables")
    .concatenate(["n_sims","drift_scaling", "threshold", "ndt", "rel_sp", "rt_multiplier"],into="inference_variables")
    .concatenate(["alt_option","self_option","PA"], into="inference_conditions")
    )
adapt_sims = adapter(simulator.sample(5))

```

adapt\_sims ends up being a dictionary with the array “inference\_conditions” of size 5 \* 3 (nbatches x n\_contextvariables), “inference\_variables” of size 5 \* 6 (nbatches x nparameters) and finally, “summary\_variables” of size 5\*1\* _2_ (nbatches x ntrials x nvariables). I am wondering whether this is right, whether having all three context variables together is “correct” and whether they should not be split into their own dimension each (similar to how the “summary\_variables” operates)? And is this even the appropriate way of incorporating context variables into the workflow in the first place, or am I completely off the mark?

1. Another concern I have is, how best to scale up the functions so that they accommodate multiple decision trials at once? Right now one batch represents a single decision trial; however, my guesstimate is that it would be more computationally efficient to simulate responses for a whole dataset (say, 96 trials) at once, with each of the 3 context variables also being generated 96 times, such that a single batch would represent a whole experiment. I can easily extend my random\_dm function to generate multiple trials - what I’m wondering is, how should I format the choice and RT variables, and especially the 3 context variables, in the adapter such that the neural network correctly interprets the relationship between a particular set of context variables (“inference conditions”) and a particular set of response variables (“summary variables”)?

Many thanks in advance!

---

<div class="post-metadata">

**Author:** ![valentin](https://yyz1.discourse-cdn.com/flex007/user_avatar/discuss.bayesflow.org/valentin/32/103_2.png) [@valentin](https://discuss.bayesflow.org/u/valentin)\
**Post date:** [August 27, 2025, 12:38pm UTC](https://discuss.bayesflow.org/t/fitting-a-value-based-decision-model-with-trial-level-context-variables/194/2 "2025-08-27T12:38:30Z")

</div>

Hi Erik,  
welcome to the forum!  
Your setting is a bit complex, so please correct me if I misread something.

Regarding your first question:  
My impression is that if you have only one combination of `alt_option`, `self_option` and `PA` for one simulation, the setup you have looks good. The `inference_conditions` basically allow you to tell the network “for this simulation, variable 1 had this value, variable 2 had this value, etc., use those as a condition without using a summary network on them”. The additional dimension is needed when you have a set of observations (i.e., multiple trials), but then you also have to process them with a summary network.

Regarding your second question:  
This is somewhat related to the first. I’m a bit uncertain about the goal you have for the simulations. For the complete dataset with 96 trials, would the participant-level variable vary between the trials, or would this stay constant? For `alt_option` and `self_option`, my guess would be that you would want to process them like the other parameters in the prior\_sampler. To simulate the whole dataset, you would then sample 96 values for each parameter in the prior, and then run the decision model for each in the `random_dm` function.  
The main difference between “inference\_conditions” and “summary\_variables” is not whether they are “context” or not, but whether they are passed through the summary network or not. If the context is the same for all trials, it makes sense to pass it as “inference\_conditions”, but if it varies you can just pass it to “summary\_variables” like the observations.

Does this help? If I misunderstood something, please let me know.

---

<div class="post-metadata">

**Author:** ![ErikStuchly](https://avatars.discourse-cdn.com/v4/letter/e/59ef9b/32.png) [@ErikStuchly](https://discuss.bayesflow.org/u/ErikStuchly)\
**Post date:** [September 8, 2025, 3:13pm UTC](https://discuss.bayesflow.org/t/fitting-a-value-based-decision-model-with-trial-level-context-variables/194/3 "2025-09-08T15:13:29Z")

</div>

Hi Valentin,

Thank you very much for your quick answer and apologies for getting back to you so late - other things have kept me busy lately.

Either way yes, to clarify, my idea was to train the model on simulated datasets of approx. 96 decision trials, since my participants should each complete around 96 trials (at least in the conditions which I care about). This means that a single call to the simulator will give me the following variables; I also describe whether they should be constant/varying within a participant (so, within a single dataset), and whether they should be constant within a batch/across batches:

- _n\_trials_ (the no. of trials of a single participant) = an integer (length=1) fixed for each participant, fixed for a batch but varying across batches
- _PA_ = a float (length=1), fixed within a participant, varying across participants within a batch/across batches
- _alt\_option_; _self\_option_ = floats (length=n\_trials), varying within a participant, varying within and across batches
- _n\_sims_; _drift\_scalings_; _thresholds_; _ndts_; _rel\_sps_; _rt\_multipliers_ = floats (length = 1), constant within a participant, varying across participants
- _accuracy_ = a 0/1 binary variable (length=n\_trials), the output variable
- _rt_ = a float (length=n\_trials), the output variable

So essentially, I create a function to generate number of trials (for all participants within a batch):

```auto
def pick_ntrials():
    return {"n_trials": random.randint(75, 96)}

```

a function to generate my decision option pairs for a single participant:

```auto
def make_optionpairs(n_trials=96):
    step = 0.05
    max_int = 19

    # [batch, n_trials]
    alt_ints = np.random.randint(0, max_int, n_trials)
    max_selfs = max_int - alt_ints - 1
    max_selfs = np.clip(max_selfs, 1, None)
    self_offsets = np.random.randint(1, max_selfs + 1) # offset at least 1
    self_ints = alt_ints + self_offsets

    alt_option = alt_ints * step
    self_option = self_ints * step

    PA = np.random.random(1) # shape: (batch_size,)

    return {
        "alt_option": alt_option, # shape [batch, n_trial]
        "self_option": self_option,
        "PA": PA
    }

```

a function to sample parameters for a single dataset/participant:

```auto
def prior_sampler():
    # Each is [batch_size]
    n_sims = np.round(np.random.uniform(1,41)) # usually one simulation per batch/participant
    drift_scaling = np.random.uniform(0.1, 10)
    threshold = np.random.uniform(0.7, 3.0)
    ndt = np.random.uniform(0,0.3)
    rel_sp = np.random.uniform(0.3,0.7)
    rt_multiplier = np.random.uniform(0.1,100)
    return {
        "n_sims": n_sims,
        "drift_scaling": drift_scaling, # shape [batch_size]
        "threshold": threshold,
        "ndt": ndt,
        "rel_sp": rel_sp,
        "rt_multiplier": rt_multiplier
    }

```

and finally, a function to generate the dataset for a single participant:

```auto
def random_dm(alt_option, self_option, PA, n_sims, drift_scaling, threshold, ndt, rel_sp, rt_multiplier, noise_constant=1, dt=0.001, max_rt=10):
    #print("self_option:", self_option, "shape:", np.shape(self_option))
    #print("PA:", PA, "shape:", np.shape(PA))
    #print("alt_option", alt_option, type(alt_option), np.shape(alt_option))
    v = (drift_scaling *
         ((1 - (alt_option - PA)**2) -
          (1 - (self_option - PA)**2))
        )
    n_sims = int(np.asarray(n_sims).item())
    #drifts = np.broadcast_to(v[:, None], (len(v), n_sims))
    drifts = np.tile(v.T, (1,n_sims))
    n_trials = drifts.shape[0]
    acc = np.full(drifts.shape, np.nan)
    rt = np.full(drifts.shape, np.nan)
    max_tsteps = int(max_rt / dt)
    
    x = np.ones(drifts.shape) * rel_sp * threshold # evidence at t=0
    ongoing = np.full(drifts.shape, True)
    
    tstep = 0
    while np.any(ongoing) and tstep < max_tsteps:
        # apply drift and noise only to ongoing trials
        moving_indices = np.where(ongoing)
        x[moving_indices] += (
            drifts[moving_indices] * dt +
            np.random.normal(loc=0, scale=noise_constant * np.sqrt(dt), size=len(moving_indices[0]))
        )
        tstep += 1
        
        # handle threshold crossing
        ended_correct = (x >= threshold)
        ended_incorrect = (x <= 0)
        
        # For correct responses
        mask = (ended_correct & ongoing)
        if np.any(mask):
            acc[mask] = 1
            rt[mask] = dt * tstep
            ongoing[mask] = False

        # For incorrect responses
        mask = (ended_incorrect & ongoing)
        if np.any(mask):
            acc[mask] = -1
            rt[mask] = dt * tstep
            ongoing[mask] = False

    rt = rt / rt_multiplier

    # Aggregate "majority" choice per trial (as in rowSums)
    row_acc = np.nansum(acc, axis=1)
    acc_averaged = np.zeros(n_trials, dtype=int)
    acc_averaged[row_acc > 0] = 1
    acc_averaged[row_acc < 0] = 0
    undecided = (row_acc == 0)
    acc_averaged[undecided] = np.random.choice([1, 0], size=np.sum(undecided))

    rt_averaged = np.nansum(rt, axis=1) + ndt

    return dict(accuracy=acc_averaged, rt=rt_averaged)

```

Then I put it together into a simulator, ensuring that pick\_ntrials is passed as a meta\_fn argument (to ensure that n\_trials stays the same within a batch, but varies across batches):

```auto
batch_size=2
simulator = bf.make_simulator(
    [make_optionpairs,prior_sampler, random_dm],
    meta_fn=pick_ntrials
)
test_sims = simulator.sample(batch_size)

```

Finally, I create an adapter like so, following your suggestion to treat the variables which are constant across trials as an _inference\_condition,_ and those which vary as _summary\_variables_:

```auto
adapter = (
    bf.Adapter()
    .to_array()
    .broadcast("n_trials", to="accuracy")
    #.constrain(["n_sims","drift_scaling", "threshold", "ndt", "rel_sp", "rt_multiplier"],lower=np.array([0,0,0,0,0,0]),upper=np.array([100,10,10,1.5,1,5]))
    .as_set(["accuracy", "rt","alt_option","self_option"])
    #.as_set(["alt_option","self_option","PA"])
    .convert_dtype("float64", "float32")
    .concatenate(["accuracy","rt","alt_option","self_option",], into="summary_variables")
    .concatenate(["n_sims","drift_scaling", "threshold", "ndt", "rel_sp", "rt_multiplier"],into="inference_variables")
    .concatenate(["PA","n_trials"], into="inference_conditions")
    )
adapt_sims = adapter(simulator.sample(3))

```

This gives me inference\_conditions of size [n\_batches,n\_vars] (where vars are _PA_ and _n\_trials_); inference\_variables of size [n\_batches,n\_vars] (where vars are the 6 parameters); and summary\_variables of size [n\_batches,n\_trials,n\_vars] (where vars are _accuracy_, _rt_, _alt\_option_, _self\_option_). To me that sounds quite right, is that the case in your opinion? 🙂

---

<div class="post-metadata">

**Author:** ![KLDivergence](https://yyz1.discourse-cdn.com/flex007/user_avatar/discuss.bayesflow.org/kldivergence/32/15_2.png) [@KLDivergence](https://discuss.bayesflow.org/u/KLDivergence)\
**Post date:** [September 11, 2025, 12:20pm UTC](https://discuss.bayesflow.org/t/fitting-a-value-based-decision-model-with-trial-level-context-variables/194/4 "2025-09-11T12:20:50Z")

</div>

Yes, the outputs shapes are reasonable!

---

<div class="post-metadata">

**Author:** ![ErikStuchly](https://avatars.discourse-cdn.com/v4/letter/e/59ef9b/32.png) [@ErikStuchly](https://discuss.bayesflow.org/u/ErikStuchly)\
**Post date:** [September 15, 2025, 9:33am UTC](https://discuss.bayesflow.org/t/fitting-a-value-based-decision-model-with-trial-level-context-variables/194/5 "2025-09-15T09:33:13Z")

</div>

Awesome, thank you very much for the confirmation!