Эксперимент Бюффона с подбрасыванием иглы: моделирование на Python

Бюффон (1707-1788)

"Опыт" Бюффона (1707—1788) с подбрасыванием иглы является, пожалуй, самым оригинальным способом оценки числа \(\pi\). В чем, собственно, его оригинальность?! Представьте, что у вас имеется игла и большой лист бумаги. Если теперь  разлиновать этот лист так, чтобы параллельные линии отстояли друг от друга на расстояние вдвое большее длины иглы, то производя всего лишь подбрасывания  иглы и отмечая факт пересечения ей линий, можно получить приближение числа \(\pi\). Неожиданно?

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

Essay of moral arithmetic
"Эссе нравственной арифметики" Бюффон, 1777

Первоначально, "эксперимент Бюффона" состоял в подбрасывании монеты, и был связан с оценкой вероятности накрытия ею различного числа ячеек при подбрасывании  над разлинованной в ячейки плоскостью. В 1812 г. П.С. Лаплас предложил использовать данный подход для оценки числа \(\pi\).

В практической реализации эксперимента с подбрасыванием иглы очень важно соблюсти "случайность" каждого броска. Ведь каждый раз подбрасывая иглу, мы в большинстве случаев держим ее определенным образом, подбрасываем приблизительно на одну и ту же высоту,  поверхность стола над которым проводится эксперимент, может иметь незначительный наклон и т.п. — все эти факторы дают свой "неслучайный" систематический вклад, в связи с чем получаемые оценки \(\pi\) получают смещение и, как результат, невысокую точность.

Опыт Бюффона с точки зрения теории вероятностей

Рассмотрим следующую задачу.

Задача Бюффона. Плоскость расчерчена параллельными прямыми, отстоящими друг от друга на расстоянии \(2h\). На плоскость наудачу брошена игла длины \(2l\). Определить вероятность того, что игла пересечет какую-либо прямую.

Определим множество элементарных событий \(\Omega\). Для этого заметим, что факт пересечения иглы с прямой характеризуется близостью центра и величиной угла между иглой и одной из параллельных прямых.

Пересечение иглы с одной из прямых наступает тогда, когда \(z<l\cdot\sin(\varphi)\), если в качестве параметров, описывающих положение иглы принять: \(\varphi\) — угол между иглой и одной из параллельных прямых; \(z\) — расстояние от центра иглы до ближайшей прямой. Учитывая, что интерес представляют только события, состоящие в пересечении иглой одной из прямых линий, множество элементарных событий \(\Omega\) будет состоять из точек прямоугольника \((\varphi,\, z)\in[0,\pi]\times[0,\, h]\). Бросание иглы в эксперименте эквивалентно случайному выбору точки из прямоугольника \(\Omega\), если  для выбранной точки с координатами \((\varphi,\, z)\) выполнено \(z<l\cdot\sin(\varphi)\), то игла пересекает прямую, в противном случае — нет. Полагая в качестве \(\mathfrak{A}\) всевозможные подмножества \(\Omega\), имеющие площадь, определим вероятность любого события \(A\in\mathfrak{A}\), как \(P(A)=\dfrac{S(A)}{S(\Omega)}\), где \(S(\cdot)\) — площадь подмножества из \(\mathfrak {A}\). Тогда вероятность события, состоящего в том, что бросок иглы закончится пересечением ее с прямой, определится выражением: $$P(z<l\cdot\sin(\varphi))=\dfrac{\int\limits_0^{\pi}l\cdot\sin
{\varphi}d\varphi}{\pi\cdot h}=\dfrac{2l}{h\pi}.
$$

Если математическая модель случайного эксперимента построена правильно, то при достаточно большом количестве испытаний отношение \(m/n\) — числа пересечений иглы с прямой к числу бросаний, будет приблизительна равна \(\dfrac{2l}{h\pi}\). Это дает возможность организовать серию экспериментов для эмпирической оценки числа \(\pi\). Например, выбрав расстояние между прямыми вдвое больше, чем длина иглы, найдем, что \(\pi\approx\dfrac{n}{m}\), где \(m\) — количество пересечений иглы хотя бы с одной из прямых линий; \(n\) — общее число подбрасываний. Для получения хороших приближений числа \(\pi\) однако необходимо проделать большое число испытаний.

Результаты отдельных экспериментов таковы: Вольф, 1850 г — 5000 испытаний (\(\pi\approx3.1596\)); Рейна, 1925 г — 2520 испытаний (\(\pi\approx3.1795\)). Архимедово приближение числа \(\pi\approx22/7\) на порядок точнее полученных оценок!

В отношении предложенного способа эмпирической оценки числа \(\pi\) правомерно поставить вопрос: сколько испытаний с подбрасыванием иглы необходимо проделать, чтобы  с высокой вероятностью можно было утверждать, что погрешность оценки числа \(\pi\)  меньше заданного порога \(\alpha\)? Такая постановка вопроса характерна для широко используемого на практике метода статистических испытаний Монте--Карло при решении сложных задач, когда единственным средством исследования явления служат результаты многократного его повторения или имитационного моделирования, например, в вычислительной среде.

Если положить \(h=2l\), то вероятность пересечения иглой прямой линии связана с искомым значением \(\pi\): \(p=1/\pi\). Зададимся точностью вычисления \(1/\pi\) — \(\alpha\) и порогом вероятности \(\varepsilon\) такими, что $$ P(\left|\dfrac{\nu_n}{n}-p\right|<\alpha)>1-\varepsilon.
$$

Задача состоит в том, чтобы определить  минимальное значение \(n\), при котором вероятность попадания частоты \(\dfrac{\nu_n}{n}\) в интервал \((p-\alpha,p+\alpha)\) больше, чем \(1-\varepsilon\). Величину \(p\), вообще говоря, следует считать неизвестной. $$\begin{array}{l} P(\left|\dfrac{\nu_n}{n}-p\right|<\alpha)= \underbrace{
\sum\limits_{n(p-\alpha)<j<n(p+\alpha)}C_n^j p^j
(1-p)^{n-j}}_{\mbox{трудно вычислить}}=\\ P(-\dfrac{\alpha\sqrt
n}{\sqrt {pq}}<\dfrac{\nu_n-np}{\sqrt{npq}}<\dfrac{\alpha\sqrt
n}{\sqrt{pq}})\approx\dfrac{1}{\sqrt{2\pi}}\int\limits_{-\alpha\dfrac{\sqrt
n}{\sqrt{pq}}}^{\alpha\dfrac{\sqrt
n}{\sqrt{pq}}}\textrm{e}^{-x^2/2}dx. \end{array}$$

Искомое число испытаний \(n\) приближенно можно найти, решив неравенство$$
\dfrac{1}{\sqrt{2\pi}}\int\limits_{-\alpha\dfrac{\sqrt
n}{\sqrt{pq}}}^{\alpha\dfrac{\sqrt
n}{\sqrt{pq}}}\textrm{e}^{-x^2/2}dx>1-\varepsilon. $$

Далее, учитывая, что в данном эксперименте значение \(\pi\) известно, и мы хотим определить \(n\) исключительно, чтобы удовлетворить свое любопытство, можно считать, что \(pq=\dfrac{1-\pi}{\pi^2}\). В реальных задачах \(p\) неизвестно (иначе зачем тогда проводить статистические испытания, если ответ итак известен). В последнем случае приходится довольствоваться лишь оценкой \(pq\leq1/4\), являющейся следствием связи \(p+q=1\).

Неравенство, представленное выше, не имеет аналитического решения в общем случае, его следует решать при конкретных значениях \(\alpha\) и \(\varepsilon\).

Пусть \(\alpha=0.001\) и \(\varepsilon=0.001\). Полагая, что \(p=1/\pi\) и численно решая неравенство, найдем \(n>232013\). Следовательно, чтобы изложенным выше способом найти оценку числа \(\pi\) с точностью до 3-го знака с вероятностью 99.9% необходимо проделать не менее 232013 подбрасываний иглы. Теперь становится ясным, почему результаты эксперимента с 5000 испытаниями привели к столь грубой оценке. Если считать \(p\) неизвестным и воспользоваться оценкой \(pq\leq1/4\), ввиду положительности подынтегральной функции, при тех же условиях, найдем \(n>267311\). Таким образом, априорное знание \(\pi\) позволяет при постановке эксперимента для достижения установленной точности сэкономить приблизительно на 35000 бросаниях иглы!

Опыт Бюффона: реализация на Python

#coding: utf-8

import numpy as  np

np.random.seed(100)

def pi_buffon(H, L, N=10**4):
    positions = np.random.uniform(0.0, H, N)
    phis = np.random.uniform(0, 2*np.pi, N)
    res = (np.abs(L*np.sin(phis)) >= positions).sum()
    return float(N)/res*2*L/H


X = np.linspace(0.1, 20, 1000) 
Y = np.linspace(0.1, 20, 1000)

xx, yy = np.meshgrid(X, Y)
zz = []
for x,y in zip(xx.ravel(), yy.ravel()): 
    zz.append(np.abs(pi_buffon(x,y) - np.pi))
    
from pylab import *

rcParams['font.family'] = 'DejaVu Sans' # Понимает русские буквы
rcParams['font.size'] = 16

contourf(xx, yy, np.log(np.array(zz).reshape(xx.shape)))
xlabel('Gap len')
ylabel('Needle len')
colorbar()
show()

В представленном фрагменте кода рассмотрена реализация эксперимента Бюффона с подбрасыванием иглы при разлиных комбинациях длины иглы и расстояний между линиями. Цветом на рисунке обозначен логарифм ошибки оценки числа \(\pi\).

Python реализация задачи Бюффона

Таким образом, наилучшая точность оценки \(\pi\) имеет место, когда длина иглы несколько меньше промежутка между линиями.
 

blog comments powered by Disqus