Вычисление интеграла методом геометрического Монте-Карло¶
Идея метода геометрического метод Монте-Карло описана в разделе geomMC.
Для реализации алгоритма зададим функцию для интегрирования:
def func(x): """ function for integration half of a circle with radius = 1 """ if abs(x) > 1: return None return (1.-x**2)**0.5
В данном случае это дуга единичной окружности с ценром в начале координат.
Предварительно зададим функцию ошибок:
def errf(prob): """ error fumction return approx. upper limit of err function for given probability """ dat="""0.00 0.00000 0.01 0.00798 0.02 0.01596 0.03 0.02393 0.04 0.03191 0.05 0.03988 3.94 0.99992 3.95 0.99992 3.96 0.99992 3.97 0.99993 3.98 0.99993 3.99 0.99993""".split('\n') dd=[[ float(j) for j in i.split()] for i in dat] for i in dd: if prob < i[1]: return i[0] return 4
Удобно взять текстовую таблицу из интернета, присвоить ее строковой переменной и автоматически превратить эти данные в функцию. Обратите внимание, что средняя часть таблицы не выведена на экран из-за экономии места. Функция возвращает верхний предел интегрирования (параметр, определяющий доверительный интервал) для заданного значения вероятности (значения интеграла).
Задаем параметры и начальные значения сумм.
# max numper of points, confidence probability and x parameter for confidence interval Nmax=100000 prob=0.98 erf=errf(prob) # number of points, sum of values and sum of squared values sum_d = 0. sum2_d = 0. nn = 0
Накапливаем число точек, сумму испытаний и сумму квадратов испытаний для быстрого вычисления среднего значения и дисперсии. Через некоторое количество испытаний обрабатываем результат и ставим точку на графике.
for i in range(Nmax): x=random() y=random() dd = 0 if y < func(x): dd = 1 sum_d += dd sum2_d += dd*dd nn += 1 if i%1000 == 0 and i > 0: av = float(sum_d)/nn disp = float(sum2_d)/nn-av**2 interv = erf*(disp/nn)**0.5 err = (av-3.14159265/4) py.errorbar(nn,err,yerr=interv) py.draw()