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

Question

Пусть 20 человек, включая ровно 3 женщины, рассаживаются случайным образом за 4 столами (обозначаются (A, B, C, D)) по 5 человек в каждой, причем все условия одинаково вероятны. Пусть X будет количеством столов, за которыми не сидят женщины. Напишите обалденное моделирование Монте-Карло, чтобы оценить ожидание X, а также оценить вероятность p , что ни одна женщина не сядет за стол A. Запустите моделирование для 3 случаев (100, 100, 10000)

Я хотел бы определить функцию, которая использует функцию random.permutation numpy для вычисления ожидаемого значения X с учетом переменного числа испытаний. Я понимаю, как это сделать на ручке и бумаге, повторяя мою коллекцию вероятностей и умножая их друг с другом, так что я могу рассчитать общую вероятность события. Это то, что я до сих пор

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)

Если я проверю свою коллекцию, это мой вывод: [0,85, 0,8421052631578947, 0,8333333333333334, 0,8235294117647058, 0,8125]

Answer 1

Реализация

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

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

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

Обсуждение

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

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

Answer 2

Вы можете умножить элементы внутри коллекции, используя 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

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

Answer 3

Вы должны явно смоделировать заседания?Если нет, то просто нарисуйте 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).