Получить число Пи из модели Эренфестов

[denny]

Хорошо известна задача Бюффона, и получение $\pi$ из серии экспериментов.
Я экспериментировал с моделью Эренфестов, и обратил внимание,
что мат ожидание времени возвращения к равновесию равно $\sqrt{\pi N/2}$, если начинать из равновесия.
То есть, мы можем получить $\pi$ и из этой модели?

Провел численные эксперименты, но результаты не стабильные, ничуть не похожи на эти.

Я не очень разбираюсь в теорвере, но базовые знания есть. Буду признателен, если на основе понятий "сходимость по вероятности", ЗБЧ и подобных помогут организовать вычисления так, чтобы из данных экспериментов на урновой модели получить $\pi$.
Или объяснить, почему это невозможно!

[svv]

denny в сообщении #1590187 писал(а):

результаты не стабильные

Вы имеете в виду, что при огромном числе испытаний точность всё равно низкая? Да, могу подтвердить. На ста миллионах испытаний результат может быть и $3.13$, и $3.15$.

Используется синтаксис C++

#include <cstdlib>
#include <ctime>
#include <iostream>
double sqr(double a) { return a*a; }
int main()
{
   using namespace std;
   srand(time(NULL));
   const int m=100000000, half=100, N=2*half;
   int left=half, equi=0;
   for (int i=0; i<m; ++i)
   {
      if (rand()%N < left) --left; else ++left;
      if (left==half) ++equi;
   }
   cout << sqr(double(m)/equi)*2/N << endl;
}

Сейчас программа выдала $3.1409$, и это ещё повезло.

Кстати, Вы знакомы с таким методом приближённого вычисления $\pi$?
Why do colliding blocks compute pi?

[Geen]

(Оффтоп)

svv в сообщении #1590216 писал(а):

программа выдала

Ужасно.... - каков период, что делать с остатком RAND_MAX по N и почему %N, а не %(N+1) ?
Не говоря уж о том, что качество rand() не гарантируется

[svv]

Geen в сообщении #1590283 писал(а):

почему %N, а не %(N+1) ?

В этой модели на каждом ходу перекладывается в другой ящик один шар, который случайно выбирается из общего множества $N$ шаров с равной вероятностью. Поэтому нужен генератор случайного номера шара ($\text{number}$) из множества $N$ целых чисел. В таком случае обычно берут диапазон $0,1,...N-1$, т.е. шары нумеруют с нуля.

Количество шаров в левом ящике хранится в переменной $\text{left}$. Если перенумеровать все шары, начиная с левого ящика, шары в левом ящике получат номера $0,...,\text{left}-1$. А условием того, что надо взять шар из левого ящика, будет $\text{number}<\text{left}.$

Замечания по поводу реализации самого ГСЧ справедливы.

[Geen]

svv в сообщении #1590315 писал(а):

В этой модели

Вы правы, я тут протупил.

[denny]

Да, спасибо, с блоками интересная модель. Но там фактически геометрия.
А мне все же непонятно, почему в вероятностной модели Бюффона есть сходимость,
причем почти как в блоках, а в урновой модели нет?
Может надо как-то по-другому среднее считать?

[sqribner48]

denny
А что же Вы хотели? Переменная equi есть случайная величина и поэтому для получения приблизительно точного значения, надо выполнить некоторое количество прогонов в каждом из которых по 1000000000 генераций случайного числа. И по результатам получить среднее значение. А Вы сделали всего один, но этого не достаточно.
Вообще то мне ваш подход не совсем был понятен. Поэтому я открыл литературу (книжку "Математика в современном мире"). Там про это было несколько по другому написано.
Есть три урны (А, В, С), в третьей урне С находятся номера шаров. В начале эксперимента все пронумерованные шары лежат в урне А. Затем последовательно, случайным образом вынимаются номера по одному. Шар с вынутым номером перекладывается из одной урны в другую, а вынутый номер возвращается обратно в свою урну. Соответственно этому описанию, я сделал алгоритм на питоне версии 3.10, который привожу ниже

Используется синтаксис Python

import random

def erenfest(m,N):
    numBalance = 0
    lstMol = [1    for i in range(0,N)]
    nAA = len(lstMol); nBB = 0
    lstNumbAA = []; lstNumbBB = []
    for i in range(0,m):
        k = random.randint(0, N-1)
        #print(k)
        if lstMol[k] == 1:
            lstMol[k] = 0
            nAA -= 1
            nBB += 1
        elif lstMol[k] == 0:
            lstMol[k] = 1
            nBB -= 1
            nAA += 1
        if nAA == nBB:
            numBalance +=  1
    return  numBalance

m = 200000000
N = 16384  # количество шаров

lstNumBal = []
pn = 10
for j in range(0,pn):
    print("--- j  ", j)
    numBalance = erenfest(m,N)
    print("numBalance  ", numBalance)
    lstNumBal.append(numBalance)
numBalSred = sum(lstNumBal)/len(lstNumBal)
print("sred numBal ", numBalSred)
pisr = (m/numBalSred)**2 *2 /N
print("pisr ", pisr)
 

В алгоритме используются следующие переменные: $nAA, nBB$ --- кол-во шаров в урне А и В соответственно;
$lstMol$ --- список меток шаров: $lstMol[k] =1$ когда k-й шар лежит в урне А и $lstMol[k] =0$ когда k-й шар лежит в урне В;
$numBalance$ --- количество случаев, когда числа шаров в урнах совпадают; $pn$---- количество выполненных экспериментов в одной серии; $m$---- число прогонов в эксперименте.
Сначала я выполнил две серии экспериментов ($m = 20000000$) по 5 в каждой, получил 3.1518 и 3.20786, затем выполнил серию из 10 экспериментов ($m = 200000000$) получил 3.1663.
Похоже, что генератор случайных чисел не очень хорош, потому что рядом стоящие числа, мало отличаются друг от друга. К сожалению повторить расчеты с большим значением $m$ или большое число серий я не могу из-за занятости да и python не позволяет. У Вас возможностей больше, поэтому как говорится, "Вам и карты в руки".

[sqribner48]

Вот, выходил из дома по делам. Чтобы не терять времени поставил программу на расчет 10-ти серий по 10 экспериментов в каждой при $N= 16834$ и $m = 20000000$. На выходе получилось: 3.11756, 3.1298, 3.1408, 3.1466, 3.1563, 3.1610, 3.1952, 3.1478, 3.1549, 3.1587. И в результате средне значение получилось 3.14407.
Конечно в других сериях это значение будет другим!
Предлагаю невозможность(?) получения числа Пи из модели Эренфестов назвать парадоксом denny-svv-Geen.
Очевидно в виду плохого качества генераторов случайных чисел, лучше использовать большое количество малых по длине серий, чем малое число очень длинных серий. (алгоритм скрибнера).
PS На сайте dxdy кто-то просил дать хорошую тему для школьных (студенческих) докладов. Мне кажется, тема поднятая denny, подойдет для таких докладов в самый раз. Там есть и теория, и программирование, и графики изменения числа шаров в урнах можно построить, и даже попробовать "повернуть историю вспять" (т.е. смоделировать возврат в исходное состояние) при малых значениях N.
И все-таки, а в самом деле, почему оно не!!??....

[Geen]

sqribner48 в сообщении #1590546 писал(а):

надо выполнить некоторое количество прогонов в каждом из которых по 1000000000 генераций случайного числа.

Плохая идея...

sqribner48 в сообщении #1590546 писал(а):

Похоже, что генератор случайных чисел не очень хорош

Вообще, мерсенновский вихрь один из лучших, в некотором классе. Вот только чем же он инициализируется?...

[sqribner48]

Geen

Geen в сообщении #1590586 писал(а):

qribner48 в сообщении #1590546
писал(а):
надо выполнить некоторое количество прогонов в каждом из которых по 1000000000 генераций случайного числа.
Плохая идея...

А чем же она плоха? (Я имел в виду не точно 1000000000 генераций, а просто некоторое количество). Разве можно по другому определить мат. ожидание случайной величины?

Geen в сообщении #1590586 писал(а):

Вообще, мерсенновский вихрь один из лучших, в некотором классе. Вот только чем же он инициализируется?...

Да я в интернете читал, что этот ГСЧ самый лучший. Насчет Ваших слов "чем же он инициализируется?" я не понял.
Если Вам не трудно, прошу разъяснить или дать ссылку. Буду бесконечно благодарен.

[Geen]

sqribner48 в сообщении #1590599 писал(а):

А чем же она плоха?

Ключевой вопрос - как инициализируется генератор.

sqribner48 в сообщении #1590599 писал(а):

прошу разъяснить

Любой ГПСЧ надо инициализировать. Обычно для этого вызывается функция seed. Если его не инициализировать явно, то что он будет выдавать... в документации на питон этого нет. А значит, зависит от версии/реализации, и в худшем случае это будет результат инициализации seed(0), то есть при последовательных запусках программы результаты будут идентичны (это бывает полезно для отладки).

sqribner48 в сообщении #1590579 писал(а):

лучше использовать большое количество малых по длине серий, чем малое число очень длинных серий

При большом периоде - нет, при инициализации "системным временем" - тоже нет (но не так категорично).

sqribner48 в сообщении #1590546 писал(а):

генератор случайных чисел не очень хорош, потому что рядом стоящие числа, мало отличаются друг от друга

При плохом инициализаторе, разве что....

[sqribner48]

Уважаемый Geen
Большое спасибо. Попробую покопаться в литературе и в интернете насчет этой seed.

[svv]

Вот ещё один подход. Принимаем ту же модель, но не используем ГПСЧ, а вычисляем вероятности на $i$-м шаге иметь определённое количество шаров в ящиках.
Пусть всего шаров $2H$. В равновесии в каждом ящике будет $H$ шаров (от слова half).
Считаем, что на $0$-м шаге у нас равновесие. На $1$-м шаге в одном из ящиков достоверно будет $H+1$ шаров, за ним далее и следим. Условимся, что элемент $p[k]$ вектора $p[]$ (где индекс меняется от $0$ до $H$) содержит вероятность того, что на текущем шаге в выбранном ящике ровно $H+k$ шаров (т.е. индекс показывает отклонение от равновесия). Тогда вероятности на следующем шаге даются матричной формулой
$\begin{bmatrix}0&\frac{H+1}{2H}&0&0&\cdots&0&0&0\\0&0&\frac{H+2}{2H}&0&\cdots&0&0&0 \\0&\frac{H-1}{2H}&0&\frac{H+3}{2H}&\cdots&0&0&0\\0&0&\frac{H-2}{2H}&0&\cdots&0&0&0 \\\vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots&\vdots \\0&0&0&0&\cdots&0&\frac{2H-1}{2H}&0\\0&0&0&0&\cdots&\frac{2}{2H}&0&1 \\0&0&0&0&\cdots&0&\frac{1}{2H}&0\end{bmatrix}\begin{bmatrix}p[0]\\p[1]\\p[2]\\p[3] \\\vdots\\p[H-2]\\p[H-1]\\p[H]\end{bmatrix}$
Левый столбец матрицы нулевой. Это значит, что вероятность равновесия $p[0]$ на текущем шаге не даёт вклад в какие-либо вероятности на следующем шаге. Система, достигшая равновесия, выбывает из игры. Ячейка $p[0]$ используется для коррекции матожидания числа шагов $A$ и сразу же очищается, чтобы учитывать только "новые поступления".

Используется синтаксис C++

#include <iostream>
int main()
{
   using namespace std;
   const int H=1000;
   double p[H+1], A=0;
   for (int k=0; k<=H; ++k)
      p[k]=0;
   p[1]=1;

   for (int i=2; i<=100000; ++i)
   {
      double q[H+1];
      for (int k=0; k<=H; ++k)
         q[k]=0;
      for (int k=1; k<=H; ++k)
         q[k-1]+=double(H+k)/(2*H)*p[k];
      for (int k=1; k<H; ++k)
         q[k+1]+=double(H-k)/(2*H)*p[k];
      for (int k=0; k<=H; ++k)
         p[k]=q[k];
      A+=p[0]*i;
   }
   double Pi=A*A/H;
   cout << Pi << endl;
}
 

Эта программа выдаёт $3.14238$.

$A$ — матожидание числа шагов до возвращения к равновесию.
$q$ — временный массив для новых значений вероятностей.
$i$ — номер шага.
$Pi$ — приближение $\pi$.

[Geen]

svv в сообщении #1590622 писал(а):

Эта программа выдаёт $3.14238$

Видимо 1000 шаров это слишком мало.

[svv]

Для $10000$ шаров $3.14167$
Для $100000$ шаров $3.1416005$
В последнем случае программа считала минут сорок, и ей потребовалось три миллиона шагов. Заранее неизвестно, сколько шагов нужно для заданной точности, поэтому лучше переделать программу так, чтобы она выдавала текущее значение, скажем, каждую секунду, ну и следить, когда оно зафиксируется. А цикл можно сделать и бесконечным.
Медленно сходится, зараза. Так хотелось получить хотя бы $3.14159$, но нет.