Solution of bulbacon's task (https://goo.gl/cCqrJw)

bulbacon_test.py

import itertools
import random
import sys
import operator
import functools
import types
import math

import numpy as np

random.seed(11)

STEPS = 5

HProb = types.SimpleNamespace(
    left=0.2,
    stay=0.4,
    right=0.4
)
HProb.array = np.array(list(vars(HProb).values()))

VProb = types.SimpleNamespace(
    up=0.3,
    stay=0.5,
    down=0.2,
)
VProb.array = np.array(list(vars(VProb).values()))

# possible directions for moving
DIRECTIONS = np.array(list(
    itertools.product([-1, 0, 1], repeat=2))
).reshape((9,2))

# probability of moving in corresponding direction
DIRECTIONS_prob = (HProb.array.reshape(-1, 1) * VProb.array).flatten()


def brute_force(iterations=1_000_000):
    print('-> i am stupid (use brute force)', end='')

    def make_step(prob_th1, prob_th2):
        rnd = random.random()
        if rnd < prob_th1:
            return -1
        elif rnd < prob_th2:
            return 0
        else:
            return 1

    count = 0

    h1, h2 = HProb.array[:1].sum(), HProb.array[:2].sum()
    v1, v2 = VProb.array[:1].sum(), VProb.array[:2].sum()

    for it in range(iterations):
        x, y = 0, 0
        if it % (iterations//10) == 0:
            print('.', end='')
            sys.stdout.flush()

        for step in range(STEPS):
            x += make_step(h1, h2)
            y += make_step(v1, v2)

        if (x, y) == (0, 0):
            count += 1

    print('')
    probability = round(count/iterations, 9)
    print('  probability', probability)
    print('')


def check_all_combinations():
    print('-> i am engineer (check all combinations)')

    probability = 0

    for path in itertools.product(range(len(DIRECTIONS)), repeat=STEPS):
        steps = DIRECTIONS[list(path)]
        steps_probability = DIRECTIONS_prob[list(path)]
        x_path, y_path = zip(*steps)

        # filter unfavourable outcomes
        if sum(x_path) != 0 or sum(y_path) != 0:
            continue

        path_probability = functools.reduce(operator.mul, steps_probability, 1)
        probability += path_probability

    print('  probability', round(probability, 9))
    print('')


def analytically():
    print('-> i am scientist (use multinomial coefficients)')

    cnt1 = math.factorial(STEPS) / (math.factorial(2) * math.factorial(2))
    cnt2 = math.factorial(STEPS) / (math.factorial(3))

    h_prob = (
        HProb.stay ** STEPS
        + cnt1 * HProb.stay * HProb.left**2 * HProb.right**2
        + cnt2 * HProb.stay**3 * HProb.left * HProb.right
    )
    v_prob = (
        VProb.stay ** STEPS
        + cnt1 * VProb.stay * VProb.up**2 * VProb.down**2
        + cnt2 * VProb.stay**3 * VProb.up * VProb.down
    )

    probability = h_prob * v_prob
    print('  probability', round(probability, 9))
    print('')


if __name__ == '__main__':
    brute_force()
    check_all_combinations()
    analytically()