Начинающий Питон Монте-Карло Симулятор - PullRequest
0 голосов
/ 28 января 2019

Я новичок в Python и работаю над упражнениями, заданными нашим инструктором.Я борюсь с этим вопросом.

В редакторе Python напишите симуляцию Монте-Карло, чтобы оценить значение числа π.В частности, выполните следующие действия: A. Создайте два массива, один с именем x, один с именем y, который содержит 100 элементов каждый, которые являются случайным и равномерно распределенными действительными числами между -1 и 1. B. Отобразите y против x в виде точек всюжет.Маркируйте свои оси соответственно.C. Запишите математическое выражение, которое определяет, какие (x, y) пары точек данных расположены в круге с радиусом 1, центрированным в (0, 0) начале графа.D. Используйте логические маски, чтобы идентифицировать точки внутри круга и наложить их на другой цвет и размер маркера поверх точек данных, которые вы уже нанесли на график в B.

Это то, что у меня есть в данный момент.

import numpy as np
import math
import matplotlib.pyplot as plt
np.random.seed(12345)
x = np.random.uniform(-1,1,100) 
y = np.random.uniform(-1,1,100) 
plt.plot(x,y) //this works


for i in x:
    newarray = (1>math.sqrt(y[i]*y[i] + x[i]*x[i]))
plt.plot(newarray)

Есть предложения?

Ответы [ 2 ]

0 голосов
/ 28 января 2019

Вы близки к решению.Я немного изменил ваш MCVE:

import numpy as np
import math
import matplotlib.pyplot as plt

np.random.seed(12345)

N = 10000
x = np.random.uniform(-1, 1, N) 
y = np.random.uniform(-1, 1, N) 

Теперь мы вычисляем критерий, который имеет смысл в этом контексте, например, расстояние от точек до начала координат:

d = x**2 + y**2

Затем мы используем Булево индексирование для различения точек внутри и за пределами Единичного круга:

q = (d <= 1)

На данный момент лежит гипотеза Монте-Карло.Мы предполагаем, что соотношение равномерно распределенных точек в круге и в плоскости U(-1,1)xU(-1,1) является репрезентативным для площади единичного круга и квадрата.Тогда мы можем статистически оценить pi = 4*(Ac/As) из соотношения точек внутри круга / квадрата.Это приводит к:

pi = 4*q.sum()/q.size # 3.1464

Наконец, мы строим график результата:

fig, axe = plt.subplots()
axe.plot(x[q], y[q], '.', color='green', label=r'$d \leq 1$')
axe.plot(x[~q], y[~q], '.', color='red', label=r'$d > 1$')
axe.set_aspect('equal')
axe.set_title(r'Monte Carlo: $\pi$ Estimation')
axe.set_xlabel('$x$')
axe.set_ylabel('$y$')
axe.legend(bbox_to_anchor=(1, 1), loc='upper left')
fig.savefig('MonteCarlo.png', dpi=120)

Это выводит:

enter image description here

0 голосов
/ 28 января 2019

, как указано в комментарии, ошибка в вашем коде: for i in x должно быть for i in xrange(len(x))

Если вы действительно хотите использовать булеву маску, как сказано в выражении, вы можете сделать что-то вроде этого

    import pandas as pd
    allpoints = pd.DataFrame({'x':x, 'y':y})

    # this is your boolean mask
    mask = pow(allpoints.x, 2) + pow(allpoints.y, 2) < 1
    circlepoints = allpoints[mask]

    plt.scatter(allpoints.x, allpoints.y)
    plt.scatter(circlepoints.x, circlepoints.y)

увеличив число точек до 10000, вы получите что-то вроде этого scatter plot

, чтобы оценить ПИ, вы можете использовать известный вывод Монте-Карло

    >>> n = 10000
    >>> ( len(circlepoints) * 4 ) / float(n)
    <<< 3.1464
...