Demonstration of a simple horseshoe map construction in response to "how would you shuffle a deck if your only source of randomness is rolling dice?" - https://old.reddit.com/r/AskStatistics/comments/1u8kgiy/an_algorithm_for_shuffling_cards_for_people_not/oseb3lb/?context=3

horseshoe_shuffle.py

import numpy as np
from scipy.stats import kendalltau
import einops
import matplotlib.pyplot as plt

def new_deck():
    return np.arange(52)
    
def roll():
#    return np.random.randint(0,26)
    d6 = np.random.randint(0,6)
    d20 = np.random.randint(0,20)
    return d6 + d20

def shuffle(deck=None, stack_order=(2,0,1)):
    if deck is None:
        deck = new_deck()
    top, rest = np.split(deck, [roll()])
    middle, bottom = np.split(rest, [roll()])
    middle = middle[::-1]
    stacks = [top, middle, bottom]
    return np.concatenate([stacks[idx] for idx in stack_order])

def simulate(n=20, stack_order=(2,0,1)):
    baseline = new_deck()
    deck = new_deck()
    scores = []
    for i in range(n):
        deck = shuffle(deck, stack_order=stack_order)
        corr = kendalltau(baseline, deck)
        scores.append(corr)
    return np.array(scores)


runs = np.array([simulate(stack_order=(2,0,1)) for _ in range(1000)])
agg1 = einops.reduce(runs, "runs sims stats -> sims stats", np.mean)

runs = np.array([simulate(stack_order=(2,1,0)) for _ in range(1000)])
agg2 = einops.reduce(runs, "runs sims stats -> sims stats", np.mean)

plt.plot(agg1)
plt.plot(agg2)
plt.legend(["statistic1","pval1","statistic2","pval2"])