Интеграция Монте-Карло с Python - PullRequest
Новая
       16

Интеграция Монте-Карло с Python

1 голос
/ 01 марта 2020

Оцените следующий интеграл с интеграцией Монте-Карло:

Я пытаюсь выполнить интеграцию Монте-Карло по приведенной ниже задаче, где p (x) - это гауссово распределение с среднее значение 1 и дисперсия 2. (см. изображение).

Мне сказали, что, как только мы берем образцы из нормального распределения, pdf исчезает в интеграле. Пожалуйста, объясните эту концепцию и как мне решить это в Python. Ниже моя попытка. 

def func(x):
return (math.exp(x))*x

mu = 1
sigma = sqrt(2)
N = 1000
areas = []
for i in range(N):
    xrand = np.zeros(N)

    for i in range (len(xrand)):
        xrand[i] = np.random.normal(mu, sigma)
        integral = 0.0

    for i in range (N):
        integral += func(xrand[i])/N

    answer = integral
    areas.append(answer)

plt.title("Distribution of areas calculated")
plt.hist(areas, 60, ec = 'black')
plt.xlabel("Areas")

integral

1 Ответ

1 голос
/ 01 марта 2020

Интеграция Монте-Карло - это способ аппроксимации сложных интегралов без вычисления их решения в замкнутой форме. Чтобы ответить на ваш вопрос, PDF исчезает, потому что все, что вам нужно сделать, это: 1) выбрать случайное значение из указанного нормального распределения, 2) вычислить значение функции в подынтегральном выражении и 3) вычислить среднее из этих значений. Обратите внимание, что PDF становится неуместным в вычислениях; это имеет значение только в том случае, если необходимо убедиться, что более вероятные значения выбираются чаще. Вы можете понимать это как взятие взвешенного среднего, если это делает вещи более интуитивными.

Вот реализация Python, основанная на вашем исходном исходном коде. 

def func(x):
    return x * math.exp(x) 

def monte_carlo(n_sample, mu, sigma):
    val_lst = []
    for _ in range(n_sample):
        x = np.random.normal(mu, sigma)
        val_lst.append(func(x))
    return mean(val_lst)

Вы можете изменить func на любую функцию по вашему выбору для выполнения аппроксимации Монте-Карло этой функции. Вы также можете редактировать параметры функции monte_carlo, если вам дано другое распределение вероятностей.

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

MAX_SAMPLE = 200 # Adjust this value as you need
x = np.arange(MAX_SAMPLE)
y = [monte_carlo(i, 1, sqrt(2)) for i in x]
plt.plot(x, y)
plt.show()

Полученный график покажет вам значение сходимости, что является приближением к значению, вычисленному из решения замкнутого вида определенного интеграла.

...