О синтаксисе вероятности Монте-Карло

Реализация

Вот наивная реализация вашей симуляции Монте-Карло. Он не рассчитан на высокую производительность, вместо этого он позволяет вам перекрестно проверить настройку и увидеть детали:

import collections
import numpy as np

def runMonteCarlo(nw=3, nh=20, nt=4, N=20):
    """
    Run Monte Carlo Simulation
    """

    def countWomen(c, nt=4):
        """
        Count Number of Women per Table
        """
        x = np.array(c).reshape(nt, -1).T  # Split permutation into tables
        return np.sum(x, axis=0)           # Sum woman per table

    # Initialization:
    comp = np.array([1]*nw + [0]*(nh-nw)) # Composition: 1=woman, 0=man
    x = []                                # Counts of tables without any woman
    p = 0                                 # Probability of there is no woman at table A  

    for k in range(N):
        c = np.random.permutation(comp)   # Random permutation, table composition
        w = countWomen(c, nt=nt)          # Count Woman per table
        nc = np.sum(w!=0)                 # Count how many tables with women 
        x.append(nt - nc)                 # Store count of tables without any woman
        p += int(w[0]==0)                 # Is table A empty?
        #if k % 100 == 0:
            #print(c, w, nc, nt-nc, p)

    # Rationalize (count->frequency)
    r = collections.Counter(x)
    r = {k:r.get(k, 0)/N for k in range(nt+1)}
    p /= N
    return r, p

Выполнение работы:

for n in [100, 1000, 10000]:
    s = runMonteCarlo(N=n)
    E = sum([k*v for k,v in s[0].items()])
    print('N=%d, P(X=k) = %s, p=%s, E[X]=%s' % (n, *s, E))

Возврат:

N=100, P(X=k) = {0: 0.0, 1: 0.43, 2: 0.54, 3: 0.03, 4: 0.0}, p=0.38, E[X]=1.6
N=1000, P(X=k) = {0: 0.0, 1: 0.428, 2: 0.543, 3: 0.029, 4: 0.0}, p=0.376, E[X]=1.601
N=10000, P(X=k) = {0: 0.0, 1: 0.442, 2: 0.5235, 3: 0.0345, 4: 0.0}, p=0.4011, E[X]=1.5924999999999998

Построение графика распределения приводит к:

import pandas as pd
axe = pd.DataFrame.from_dict(s[0], orient='index').plot(kind='bar')
axe.set_title("Monte Carlo Simulation")
axe.set_xlabel('Random Variable, $X$')
axe.set_ylabel('Frequency, $F(X=k)$')
axe.grid()

Расхождение с альтернативной версией

Внимание: этот метод не решает поставленную задачу!

Если мы реализуем другую версию моделирования, в которой мы изменим способ проведения случайного эксперимента следующим образом:

import random
import collections

def runMonteCarlo2(nw=3, nh=20, nt=4, N=20):
    """
    Run Monte Carlo Simulation
    """

    def one_experiment(nt, nw):
        """
        Table setup (suggested by @Inon Peled)
        """
        return set(random.randint(0, nt-1) for _ in range(nw)) # Sample nw times from 0 <= k <= nt-1

    c = collections.Counter()             # Empty Table counter
    p = 0                                 # Probability of there is no woman at table A  

    for k in range(N):
        exp = one_experiment(nt, nw)      # Select table with at least one woman
        c.update([nt - len(exp)])         # Update Counter X distribution
        p += int(0 not in exp)            # There is no woman at table A (table 0)

    # Rationalize:
    r = {k:c.get(k, 0)/N for k in range(nt+1)}
    p /= N

    return r, p

Он возвращает:

N=100, P(X=k) = {0: 0.0, 1: 0.41, 2: 0.51, 3: 0.08, 4: 0.0}, p=0.4, E[X]=1.67
N=1000, P(X=k) = {0: 0.0, 1: 0.366, 2: 0.577, 3: 0.057, 4: 0.0}, p=0.426, E[X]=1.691
N=1000000, P(X=k) = {0: 0.0, 1: 0.37462, 2: 0.562787, 3: 0.062593, 4: 0.0}, p=0.42231, E[X]=1.687973

Эта вторая версия сходится к другим значениям, и она явно не эквивалентна первой версии, она не отвечает на тот же вопрос.

Обсуждение

Чтобы определить, какая реализация является правильной, я вычислил выборочные пространства и вероятности для обеих реализаций. Кажется, первая версия верна, потому что она учитывает, что вероятность того, что женщина сядет за стол, зависит от того, кто был выбран ранее. Во второй версии это не учитывается, поэтому не нужно знать, сколько людей и сколько человек может сесть за стол.

Это хорошая проблема, потому что оба ответа дают близкие результаты. Важной частью работы является правильная настройка входов Монте-Карло.

Вы можете умножать элементы внутри коллекции, используя functools.reduce в Python 3.x.

from functools import reduce
event_probability = reduce(lambda x, y: x*y, collection)

Итак, в вашем коде:

from functools import reduce

T = 4       # number of tables
N = 20      # number of persons. Assumption: N is a multiple of T.
K = 5       # capacity per table
W = 3       # number of women. Assumption: first W of N persons are women.
M = 100      #number of trials

collection = []

for i in range(K):
    x = (((N-W)-i)/(N-i))
    collection.append(x)

event_probability = reduce(lambda x, y: x*y, collection)

print(collection)
print(event_probability)

Выход:

[0.85, 0.8421052631578947, 0.8333333333333334, 0.8235294117647058, 0.8125] # collection
0.3991228070175438 # event_probability

Затем вы можете использовать результат для завершения своего кода.

Вам нужно явно моделировать сеансы? Если нет, то просто нарисуйте 3 раза наугад с заменой от 1..4, чтобы имитировать одно сидение, то есть:

def one_experiment():
    return set(random.randint(1, 4) for _ in range(3))  # Distinct tables with women.

Затем желаемые значения получаются следующим образом, где N - количество экспериментов для любого случая.

expectation_of_X = sum(4 - len(one_experiment()) for _ in range(N)) / float(N)
probability_no_women_table_1 = sum(1 not in one_experiment() for _ in range(N)) / float(N)

Для большого N получаемые значения должны быть примерно p = (3/4) ^ 3 и E [X] = (3 ^ 3) / (4 ^ 2).