Code to reproduce the simulation of André and de Langhe's response to Walasek, Mullett and Stewart.

Bets_S1.csv

simulate_choices.py

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns

# Reading the bets from S1 of W&S: 64 bets for each condition
bets_s1 = pd.read_csv("Bets_S1.csv")

# Unique ID for each bet
bets_s1["bet"] = bets_s1[["gain", "loss"]].apply(
    lambda x: f"{x.gain:.0f}-{x.loss:.0f}", axis=1
)
# Computing the expected value of each bet
bets_s1["EV"] = bets_s1["gain"] - bets_s1["loss"]


def simulate_choices(df, N=1000):
    """
    Simulate the choices for N participants: They see all 64 bets in a random order,
    and accept 50% of the bets, favoring those with the highest expected value.
    """
    all_participants = []
    for pid in range(N):  # Simulate one participant
        participant = df[["bet", "gain", "loss", "EV"]].sample(
            n=64, replace=False
        )  # We sample all the gambles, without replacement, in a random order
        choices = []
        for i in range(64):  # For each of the gamble
            prev_evs = participant.iloc[0:i][
                "EV"
            ]  # Expected value of all the previously encountered gambles
            current_ev = participant.iloc[i]["EV"]  # Expected value of the focal gamble
            ap = (
                current_ev > prev_evs
            ).mean()  # Acceptance probability: Proportion of previously encountered gambles that have a smaller EV than the focal gamble
            ap = (
                1 if np.isnan(ap) else ap
            )  # If it is the first gamble, it is accepted with probability 1.
            rp = 1 - ap  # Rejection probability
            choices.append(
                np.random.choice([0, 1], p=[rp, ap])
            )  # The choice is recorded
        participant["choices"] = choices
        participant["pid"] = pid
        all_participants.append(participant)
    return pd.concat(all_participants)  # Return the choices of all participants


np.random.seed(189123789)  # Seed for reproducibility

# Simulate 1000 choices in each condition
t0 = simulate_choices(bets_s1[bets_s1.condition == 0], 1000)
t0["condition"] = 0
t1 = simulate_choices(bets_s1[bets_s1.condition == 1], 1000)
t1["condition"] = 1
t2 = simulate_choices(bets_s1[bets_s1.condition == 2], 1000)
t2["condition"] = 2
t3 = simulate_choices(bets_s1[bets_s1.condition == 3], 1000)
t3["condition"] = 3

# All choices
choices = pd.concat([t0, t1, t2, t3]).reset_index(drop=True)

# Choices for all common gambles
choices_common = choices[
    choices.gain.apply(lambda x: x in [12, 16, 20])
    & choices.loss.apply(lambda x: x in [12, 16, 20])
]

# Plotting the choices
sns.set_context("paper")
g = sns.FacetGrid(
    col="condition", col_wrap=2, data=t, height=6.38 / 2, sharex=True, sharey=True
)
axes = g.axes.flatten()
for i, ax in enumerate(axes):
    # Summary of acceptance probabilities for all gain/loss pairs
    data = (
        choices_common[choices_common.condition == i]
        .groupby(["gain", "loss"])
        .mean()
        .reset_index()
        .pivot("gain", "loss", "choices")
    )
    # Heatmap of acceptance probabilities
    sns.heatmap(
        data, ax=ax, cbar=False, vmin=0, vmax=1, cmap="RdBu", annot=True, fmt=".0%"
    )
    ax.invert_yaxis()
    ax.set_xticklabels([12, 16, 20])
    ax.set_yticklabels([12, 16, 20])
    ax.tick_params(axis="both", which="both", length=0)
    
# Ticks and labels
axes[0].set_xlabel("")
axes[0].set_ylabel("Prospective Gain")
axes[0].set_title("G20-L20")
axes[1].set_xlabel("")
axes[1].set_ylabel("")
axes[1].set_title("G20-L40")
axes[2].set_ylabel("Prospective Gain")
axes[2].set_xlabel("Prospective Loss")
axes[2].set_title("G40-L20")
axes[3].set_ylabel("")
axes[3].set_xlabel("Prospective Loss")
axes[3].set_title("G40-L40")
plt.tight_layout()
plt.savefig("SimulatedAcceptanceRate.png", dpi=100)