def sim_bios(n, doubleNVRAM=False, D=45):
    deck = [1]*n + [0]*(D-n)
    #NVRAM draw
    draw = random.sample(deck, 6)
    if 1 in draw:
        return 1
    D -= 4
    #first draw
    deck = [1]*n + [0]*(D-n)
    draw = random.sample(deck, 5)
    if 1 in draw:
        return 1

    if doubleNVRAM:
        draw = random.sample(deck, 6)
        if 1 in draw:
            return 1
        D -= 4

    #second draw
    deck = [1]*n + [0]*(D-n)
    draw = random.sample(deck, 5)
    if 1 in draw:
        return 1
    return 0

def sim_andy(n, mulligan=True, D=45):
    if sim_draw(3, 45, 9):
        return 1
    if mulligan:
        return sim_draw(3, 45, 9) == 1
    else:
        return 0



andy = 0.0
andyNoMull = 0.0
bios = 0.0
doubleBios = 0.0
trials = 0
for i in range(1000000):
    andy += sim_andy(3)
    andyNoMull += sim_andy(3, False)
    bios += sim_bios(3, False)
    doubleBios += sim_bios(3, True)
    trials += 1

print "andy %0.4f andyNoMull %0.4f bios %0.4f doubleBios %0.4f"%(
  andy/trials,
  andyNoMull/trials,
  bios/trials,
  doubleBios/trials)

> $ python  tmp.py
andy 0.7463 andyNoMull 0.4975 bios 0.7108 doubleBios 0.8310