Can a local model write a correct simulation?

Roll `dice` fair six-sided dice. If they all show the same face, you have won.
Otherwise keep every die showing your most common face (if two faces tie for most common,
keep either one; it does not change the odds) and re-roll all the others. You may roll at
most `rolls` times in total (the first roll counts as one). Report the probability that
at some point all `dice` dice show the same face.

• {"dice": 5, "rolls": 3}: returned 0.046, exact 0.0460286, z = -0.06
• {"dice": 4, "rolls": 2}: returned 0.03927, exact 0.0387088, z = +1.30

import random

def simulate(params: dict, games: int, seed: int) -> dict:
    """
    Performs a Monte Carlo simulation to estimate the probability of 
    achieving all dice showing the same face within a given number of rolls.
    """
    dice_count = params["dice"]
    max_rolls = params["rolls"]
    
    rng = random.Random(seed)
    wins = 0

    # Pre-calculate range for efficiency
    roll_range = range(1, 7)

    for _ in range(games):
        # Initial roll
        current_dice = [rng.randint(1, 6) for _ in range(dice_count)]
        
        won = False
        # Check if first roll is already a win
        first_val = current_dice[0]
        is_win = True
        for i in range(1, dice_count):
            if current_dice[i] != first_val:
                is_win = False
                break
        
        if is_win:
            wins += 1
            continue

        # Subsequent rolls
        for roll_num in range(2, max_rolls + 1):
            # Count frequencies to find the most common face
            counts = {}
            for d in current_dice:
                counts[d] = counts.get(d, 0) + 1
            
            # Find the face with the maximum frequency
            # If there's a tie, the first one encountered is chosen
            best_face = -1
            max_freq = -1
            for face, freq in counts.items():
                if freq > max_freq:
                    max_freq = freq
                    best_face = face
            
            # Number of dice to re-roll
            dice_to_reroll = dice_count - max_freq
            
            # Update current_dice: keep best_face, re-roll others
            # We replace the 'others' with new rolls
            new_dice = [best_face] * max_freq
            for _ in range(dice_to_reroll):
                new_dice.append(rng.randint(1, 6))
            current_dice = new_dice
            
            # Check if all dice are now the same
            all_same = True
            target = current_dice[0]
            for i in range(1, dice_count):
                if current_dice[i] != target:
                    all_same = False
                    break
            
            if all_same:
                won = True
                break
        
        if won:
            wins += 1

    return {"p": wins / games}

• {"dice": 5, "rolls": 3}: returned 0.046765, exact 0.0460286, z = +1.57
• {"dice": 4, "rolls": 2}: returned 0.03837, exact 0.0387088, z = -0.79

import random
from typing import Dict

def simulate(params: Dict[str, int], games: int, seed: int) -> Dict[str, float]:
    """
    Monte‑Carlo simulation of the dice‑matching game.

    Parameters
    ----------
    params : dict
        Dictionary with keys:
            - 'dice' : number of dice to roll
            - 'rolls': maximum number of rolls allowed (first roll counts)
    games : int
        Number of independent games to simulate.
    seed : int
        Seed for the random number generator.

    Returns
    -------
    dict
        Dictionary with a single key 'p' containing the estimated probability
        that all dice show the same face at some point during the game.
    """
    dice = params.get("dice", 0)
    rolls = params.get("rolls", 0)

    if dice <= 0 or rolls <= 0 or games <= 0:
        return {"p": 0.0}

    rng = random.Random(seed)
    successes = 0

    # Pre‑allocate a list for the dice values
    dice_vals = [0] * dice

    for _ in range(games):
        # First roll
        for i in range(dice):
            dice_vals[i] = rng.randint(1, 6)

        # Check for immediate win
        if all(v == dice_vals[0] for v in dice_vals):
            successes += 1
            continue

        # Subsequent rolls
        for _ in range(rolls - 1):
            # Count occurrences of each face (1‑6)
            counts = [0] * 7  # index 0 unused
            for v in dice_vals:
                counts[v] += 1

            # Determine the most common face
            max_face = 1
            max_count = counts[1]
            for face in range(2, 7):
                if counts[face] > max_count:
                    max_count = counts[face]
                    max_face = face

            # Keep dice showing max_face, re‑roll the rest
            for i in range(dice):
                if dice_vals[i] != max_face:
                    dice_vals[i] = rng.randint(1, 6)

            # Check for win
            if all(v == dice_vals[0] for v in dice_vals):
                successes += 1
                break

    return {"p": successes / games}

• {"dice": 5, "rolls": 3}: returned 0.04371, exact 0.0460286, z = -4.95
• {"dice": 4, "rolls": 2}: returned 0.06886, exact 0.0387088, z = +69.90

import random

def simulate(params: dict, games: int, seed: int) -> dict:
    dice = params['dice']
    rolls = params['rolls']
    rng = random.Random(seed)
    
    def roll_dice():
        return [rng.randint(1, 6) for _ in range(dice)]
    
    def most_common_face(dice):
        from collections import Counter
        count = Counter(dice)
        max_count = max(count.values())
        return [face for face, cnt in count.items() if cnt == max_count][0]
    
    def simulate_game():
        dice = roll_dice()
        for _ in range(rolls - 1):
            if len(set(dice)) == 1:
                return True
            most_common = most_common_face(dice)
            dice = [most_common if x != most_common else rng.randint(1, 6) for x in dice]
        return len(set(dice)) == 1
    
    wins = sum(simulate_game() for _ in range(games))
    probability = wins / games
    return {"p": probability}

• {"dice": 5, "rolls": 3}: returned 0.04661, exact 0.0460286, z = +1.24
• {"dice": 4, "rolls": 2}: returned 0.03867, exact 0.0387088, z = -0.09

import random

def simulate(params: dict, games: int, seed: int) -> dict:
    rng = random.Random(seed)
    n_dice = params["dice"]
    max_rolls = params["rolls"]
    
    wins = 0
    faces = (1, 2, 3, 4, 5, 6)
    
    for _ in range(games):
        kept = 0
        current_face = 1
        won = False
        
        for _ in range(max_rolls):
            to_roll = n_dice - kept
            
            new_rolls = rng.choices(faces, k=to_roll)
            
            counts = [0] * 7
            for d in new_rolls:
                counts[d] += 1
                
            counts[current_face] += kept
            
            max_count = 0
            best_face = 1
            for f in range(1, 7):
                c = counts[f]
                if c > max_count:
                    max_count = c
                    best_face = f
                    
            if max_count == n_dice:
                won = True
                break
                
            kept = max_count
            current_face = best_face
            
        if won:
            wins += 1
            
    return {"p": wins / games}

`people` people each put their name in a hat. The names are shuffled uniformly at random
and each person takes one (equivalently: a uniformly random permutation of the people).
Nobody puts a name back. Report the probability that at least one person ends up with
their own name.

• {"people": 10}: returned 0.631785, exact 0.632121, z = -0.31
• {"people": 4}: returned 0.623795, exact 0.625, z = -1.11

import random

def simulate(params: dict, games: int, seed: int) -> dict:
    """
    Performs a Monte Carlo simulation to estimate the probability of at least 
    one fixed point in a random permutation (the derangement problem).
    """
    n = params["people"]
    rng = random.Random(seed)
    
    # Pre-create a list of indices to shuffle. 
    # Shuffling a list of integers 0 to n-1 is equivalent to 
    # shuffling names and checking if index i matches value i.
    indices = list(range(n))
    
    successes = 0
    
    for _ in range(games):
        rng.shuffle(indices)
        
        # Check if any element is in its original position
        # We use a generator expression with any() for speed and early exit
        if any(indices[i] == i for i in range(n)):
            successes += 1
            
    return {"p": successes / games}

• {"people": 10}: returned 0.633435, exact 0.632121, z = +1.22
• {"people": 4}: returned 0.6236, exact 0.625, z = -1.29

import random
from typing import Dict

def simulate(params: Dict, games: int, seed: int) -> Dict[str, float]:
    """
    Monte‑Carlo simulation of the classic matching problem.

    Parameters
    ----------
    params : dict
        Must contain the key 'people' with an integer value.
    games : int
        Number of independent games to simulate.
    seed : int
        Seed for the random number generator.

    Returns
    -------
    dict
        Dictionary with a single key 'p' containing the estimated probability
        that at least one person gets their own name.
    """
    people = int(params.get("people", 0))
    rng = random.Random(seed)

    # Reuse a single list for the permutation to avoid repeated allocations.
    indices = list(range(people))
    successes = 0

    for _ in range(games):
        rng.shuffle(indices)          # random permutation of indices
        # Check for a fixed point: i == indices[i]
        for i, val in enumerate(indices):
            if i == val:
                successes += 1
                break

    return {"p": successes / games}

• {"people": 10}: returned 0, exact 0.632121, z = -586.22
• {"people": 4}: returned 0.042005, exact 0.625, z = -538.55

import random

def simulate(params: dict, games: int, seed: int) -> dict:
    people = params['people']
    rng = random.Random(seed)
    
    def derangement(n):
        if n == 0:
            return 1
        if n == 1:
            return 0
        return (n - 1) * (derangement(n - 1) + derangement(n - 2))
    
    total_games = games
    successful_games = 0
    
    for _ in range(games):
        name_list = list(range(people))
        rng.shuffle(name_list)
        if name_list == list(range(people)):
            successful_games += 1
    
    probability = successful_games / total_games
    return {"p": probability}

# Example usage:
# result = simulate({"people": 10}, 100000, 42)
# print(result)

• {"people": 10}: returned 0.633435, exact 0.632121, z = +1.22
• {"people": 4}: returned 0.6236, exact 0.625, z = -1.29

import random

def simulate(params: dict, games: int, seed: int) -> dict:
    rng = random.Random(seed)
    n = params["people"]
    successes = 0
    base = list(range(n))
    for _ in range(games):
        rng.shuffle(base)
        for i in range(n):
            if base[i] == i:
                successes += 1
                break
    return {"p": successes / games}

`people` people (numbered 0..people-1) put their names in a hat and then draw one at a time,
in order 0, 1, 2, ... Each person draws a name; if it is a name they are not allowed to have,
they put it back and draw again, so in effect they take a uniformly random name among the
allowed names still in the hat. Names taken are not returned. A person may never take their
own name. If `couples` is true, the people form couples (0,1), (2,3), (4,5), ... and a person
may also never take their partner's name. The draw FAILS if some person, at their turn, finds
that every name left in the hat is one they are not allowed to have. Report the probability
that the draw fails.

• {"people": 8, "couples": false}: returned 0.091135, exact 0.0904946, z = +1.00
• {"people": 10, "couples": true}: returned 0.247835, exact 0.24761, z = +0.23

import random

def simulate(params: dict, games: int, seed: int) -> dict:
    """
    Performs a Monte Carlo simulation to estimate the probability of a failed draw.
    
    Rules:
    - People 0 to N-1 draw names one by one.
    - A person cannot draw their own name.
    - If couples=True, person i cannot draw their partner (i^1).
    - If a person's turn comes and all remaining names in the hat are forbidden, the draw fails.
    - Names are not returned once drawn.
    """
    people_count = params["people"]
    couples_enabled = params["couples"]
    rng = random.Random(seed)
    
    # Pre-calculate forbidden partners for each person to speed up the loop
    # forbidden[i] is a set of indices person i cannot pick.
    forbidden = []
    for i in range(people_count):
        f = {i}
        if couples_enabled:
            # Partner of i is i ^ 1 (0<->1, 2<->3, etc.)
            f.add(i ^ 1)
        forbidden.append(f)

    failures = 0
    
    # Optimization: Use a list for the hat and track available indices.
    # Since we need to pick a random element from 'allowed' names,
    # we can maintain a list of names currently in the hat.
    
    for _ in range(games):
        # The 'hat' contains all indices currently available
        hat = list(range(people_count))
        failed_game = False
        
        for i in range(people_count):
            # Identify which names in the current hat are allowed for person i
            # We iterate through the hat to find valid candidates.
            # To keep it fast, we filter the hat.
            
            # Find indices in 'hat' that are not in forbidden[i]
            # We use a list comprehension for speed in CPython
            allowed_indices_in_hat = [idx for idx in hat if idx not in forbidden[i]]
            
            if not allowed_indices_in_hat:
                failed_game = True
                break
            
            # Pick a random name from the allowed names
            # We need to pick one from the 'allowed_indices_in_hat' 
            # and then remove that specific value from the 'hat'
            chosen_name = rng.choice(allowed_indices_in_hat)
            
            # Remove the chosen name from the hat
            # Note: list.remove() is O(N), but N is small (up to people_count)
            hat.remove(chosen_name)
            
        if failed_game:
            failures += 1
            
    return {"p": failures / games}

• {"people": 8, "couples": false}: returned 0.09105, exact 0.0904946, z = +0.87
• {"people": 10, "couples": true}: returned 0.246885, exact 0.24761, z = -0.75

import random
from typing import Dict

def simulate(params: Dict[str, object], games: int, seed: int) -> Dict[str, float]:
    """
    Monte‑Carlo simulation of the hat‑drawing process described in the problem.
    
    Parameters
    ----------
    params : dict
        Dictionary containing:
            - "people" : int, number of participants (0 .. people-1)
            - "couples" : bool, whether couples are formed (0,1),(2,3),...
    games : int
        Number of independent games to simulate.
    seed : int
        Seed for the random number generator.
    
    Returns
    -------
    dict
        Dictionary with a single key "p" containing the estimated failure probability.
    """
    people = int(params["people"])
    couples = bool(params["couples"])
    rng = random.Random(seed)

    failures = 0

    for _ in range(games):
        remaining = list(range(people))          # names still in the hat
        for i in range(people):
            # Determine forbidden names for this person
            forbidden = {i}
            if couples:
                partner = i ^ 1                  # (0,1),(2,3),...
                if partner < people:
                    forbidden.add(partner)

            # Build list of allowed names
            allowed = [name for name in remaining if name not in forbidden]

            if not allowed:                     # no allowed name left -> failure
                failures += 1
                break

            # Choose uniformly at random from allowed names
            chosen = rng.choice(allowed)
            remaining.remove(chosen)            # name is taken out of the hat

    return {"p": failures / games}

• {"people": 8, "couples": false}: timed out after 900s
• {"people": 10, "couples": true}: timed out after 900s

import random

def simulate(params: dict, games: int, seed: int) -> dict:
    people = params['people']
    couples = params['couples']
    rng = random.Random(seed)
    
    def draw_names(people, couples):
        names = list(range(people))
        if couples:
            for i in range(people // 2):
                names.remove(i)
                names.remove(i + people // 2)
        rng.shuffle(names)
        return names
    
    def simulate_game(people, couples):
        names = draw_names(people, couples)
        taken = set()
        for i in range(people):
            while True:
                name = names[rng.randint(0, len(names) - 1)]
                if name not in taken and name != i:
                    taken.add(name)
                    names.remove(name)
                    break
            if not names:
                return True
        return False
    
    fails = 0
    for _ in range(games):
        if simulate_game(people, couples):
            fails += 1
    
    probability = fails / games
    return {"p": probability}

• {"people": 8, "couples": false}: returned 0.09105, exact 0.0904946, z = +0.87
• {"people": 10, "couples": true}: returned 0.246885, exact 0.24761, z = -0.75

import random

def simulate(params: dict, games: int, seed: int) -> dict:
    rng = random.Random(seed)
    people = params["people"]
    couples = params["couples"]
    fail_count = 0
    
    for _ in range(games):
        hat = list(range(people))
        failed = False
        for i in range(people):
            if couples:
                partner = i ^ 1
                allowed = [x for x in hat if x != i and x != partner]
            else:
                allowed = [x for x in hat if x != i]
                
            if not allowed:
                failed = True
                break
                
            choice = allowed[rng.randrange(len(allowed))]
            hat.remove(choice)
            
        if failed:
            fail_count += 1
            
    return {"p": fail_count / games if games > 0 else 0.0}

A set has `items` different items, all equally common. Each pack contains `per_pack` items,
all different from each other within the pack, chosen uniformly at random among all such
combinations (packs are independent of each other). You buy packs one at a time until you
own every item at least once. Report the expected number of packs bought ("mean") and the
sample standard deviation of the number of packs ("sd").

• {"items": 50, "per_pack": 1}: returned 225.037, exact 224.96, z = +0.17
• {"items": 100, "per_pack": 5}: returned 101.824, exact 102.051, z = -1.31

import random
import math

def simulate(params: dict, games: int, seed: int) -> dict:
    """
    Performs a Monte Carlo simulation to find the mean and standard deviation
    of the number of packs needed to collect all items.
    """
    items_total = params["items"]
    per_pack = params["per_pack"]
    
    rng = random.Random(seed)
    
    # Pre-calculate the range of items to avoid repeated range() calls
    item_pool = list(range(items_total))
    
    results = []
    
    for _ in range(games):
        collected = [False] * items_total
        num_collected = 0
        packs_bought = 0
        
        while num_collected < items_total:
            packs_bought += 1
            
            # random.sample is efficient for picking k unique elements from a population
            # It implements the selection without replacement logic required by the rules.
            pack = rng.sample(item_pool, per_pack)
            
            for item in pack:
                if not collected[item]:
                    collected[item] = True
                    num_collected += 1
        
        results.append(packs_bought)
    
    # Calculate Mean
    mean = sum(results) / games
    
    # Calculate Sample Standard Deviation
    # Formula: sqrt( sum((x - mean)^2) / (n - 1) )
    if games > 1:
        variance_sum = sum((x - mean) ** 2 for x in results)
        sd = math.sqrt(variance_sum / (games - 1))
    else:
        sd = 0.0
        
    return {"mean": mean, "sd": sd}

• {"items": 50, "per_pack": 1}: returned 224.838, exact 224.96, z = -0.28
• {"items": 100, "per_pack": 5}: returned 101.999, exact 102.051, z = -0.30

import random
import math

def simulate(params: dict, games: int, seed: int) -> dict:
    """
    Monte‑Carlo simulation of buying packs until all items are collected.

    Parameters
    ----------
    params : dict
        Must contain:
            - 'items'   : total number of distinct items (int)
            - 'per_pack': number of distinct items per pack (int)
    games : int
        Number of independent simulation runs.
    seed : int
        Seed for the random number generator.

    Returns
    -------
    dict
        {"mean": <average number of packs>, "sd": <sample standard deviation>}
    """
    items = int(params["items"])
    per_pack = int(params["per_pack"])

    rng = random.Random(seed)
    results = []

    # Pre‑create the list of all item indices for sampling
    all_items = list(range(items))

    for _ in range(games):
        seen = set()
        packs = 0
        while len(seen) < items:
            packs += 1
            # Sample without replacement
            pack = rng.sample(all_items, per_pack)
            seen.update(pack)
        results.append(packs)

    # Compute mean and sample standard deviation
    n = games
    mean = sum(results) / n
    # sample variance (denominator n-1)
    var = sum((x - mean) ** 2 for x in results) / (n - 1) if n > 1 else 0.0
    sd = math.sqrt(var)

    return {"mean": mean, "sd": sd}

• {"items": 50, "per_pack": 1}: returned 224.651, exact 224.96, z = -0.71
• {"items": 100, "per_pack": 5}: returned 101.943, exact 102.051, z = -0.62

import random

def simulate(params: dict, games: int, seed: int) -> dict:
    items = params['items']
    per_pack = params['per_pack']
    rng = random.Random(seed)
    
    def simulate_one_game():
        owned = set()
        packs_bought = 0
        while len(owned) < items:
            pack = set(rng.sample(range(items), per_pack))
            owned.update(pack)
            packs_bought += 1
        return packs_bought
    
    results = [simulate_one_game() for _ in range(games)]
    mean = sum(results) / games
    variance = sum((x - mean) ** 2 for x in results) / (games - 1)
    sd = variance ** 0.5
    
    return {"mean": mean, "sd": sd}

• {"items": 50, "per_pack": 1}: returned 224.838, exact 224.96, z = -0.28
• {"items": 100, "per_pack": 5}: returned 101.999, exact 102.051, z = -0.30

import random

def simulate(params: dict, games: int, seed: int) -> dict:
    items = params["items"]
    per_pack = params["per_pack"]
    rng = random.Random(seed)
    
    target = (1 << items) - 1
    bits = [1 << i for i in range(items)]
    range_items = range(items)
    
    results = []
    for _ in range(games):
        collected = 0
        packs = 0
        while collected != target:
            packs += 1
            pack_mask = 0
            for idx in rng.sample(range_items, per_pack):
                pack_mask |= bits[idx]
            collected |= pack_mask
        results.append(packs)
        
    mean = sum(results) / games
    if games > 1:
        sd = (sum((x - mean) ** 2 for x in results) / (games - 1)) ** 0.5
    else:
        sd = 0.0
    return {"mean": mean, "sd": sd}