Управление рисками на примере Санкт-Петербургского парадокса

[sic]

UPD: игнорируйте это, я решал немного другую задачу!

Скрытый текст

Вот как раз после исходной статьи стало интересно провести небольшое моделирование, и я удивлен, что у Вас получилось B_mean = 16.52, а у меня в среднем 36 (на 2^32 исходов). Даже для ста миллионов игр, уже в среднем 30 получается. Возможно ошибка в:

B[offset:offset + size] = -np.log2(rng.random(size))

log2 не находит положение первого взведенного бита корректно, ну для примера log2(1000b) очень близок к log2(0111b) но вот только в первой игре выигрышь все 32 рубля, а во второй всего 16 (или 0, если считать trailing нули, а не ведущие). На самом деле я не вполне понимаю еще из-за наличия минуса перед логарифмом, но в любом случае 2 ** B не соответствует выигрышу никаким образом. Ожидается что-то похожее на:

static long PlayGame()
{
    ulong x = rng.NextUInt64();
    int tosses = TrailingZerosBinary(x);
    return 1L << (tosses + 1); // prize value
}

Скрытый текст

Совсем говнокод, у меня древний .NET и нет нормальной TrailingZeros, но как-то так:

static int TrailingZerosBinary(ulong x)
    {
        if (x == 0) return 64;
        int n = 63;
        if ((x & 0x00000000FFFFFFFFUL) > 0) { n -= 32; } else x >>= 32;
        if ((x & 0x000000000000FFFFUL) > 0) { n -= 16; } else x >>= 16;
        if ((x & 0x00000000000000FFUL) > 0) { n -= 8; } else x >>= 8;
        if ((x & 0x000000000000000FUL) > 0) { n -= 4; } else x >>= 4;
        if ((x & 0x0000000000000003UL) > 0) { n -= 2; } else x >>= 2;
        return n - (int)(x & 1);
    }

Еще раз все у себя перепроверил, в среднем B_meanуже после 20000 игр достигает вашего значения (но на таком небольшом количестве игр оно не очень хорошо стабилизируется). По сути конечно мало меняет, а по факту - много.

[celen]

Не может быть у вас средний выигрыш 36 для серии из 2^32 исходов. Это невозможный уровень удачи. Может быть, у вас другая модель, например вы "выплачиваете" в своей симуляции агенту вдвое больший выигрыш, чем я? Тогда это то же что и выигрыш у меня 18, что в принципе допустимо.

На серии в 20000 ~ 2^14 по моей модели ожидается средний выигрыш 7 +- 3. Выход за этот диапазон явно показывает численную несогласованность наших моделей. Но если вы достигли моего значения (16) выплачивая игроку вдвое больше, то на самом деле вы достигли 8.

[sic]

Это так, я взял игру из статьи на которую Вы ссылаетесь, а там стартовый выигрыш 2 рубля :) Гоняю теперь с таким же условием, что у Вас, теперь в среднем что-то близкое получается, можно сказать тоже самое, потому что от запуска к запуску в диапазоне 10% гуляет. Теоретически на некоторых запусках совсем все что угодно может быть.

Скрытый текст

Буквально на десятый прогон словил джекпот.

[10^10] Игр: 4 100 000 000
Средний выигрыш: $153,02
Максимальный выигрыш: $549 755 813 888

Если этот выброс исключить (вот насколько реален кстати), то да, средний средний выигрыш это 18 с небольшим.

Но я не понял этого трюка с логарифмом все равно. Он же вещественный, как-то должен к целому приводиться, отбрасыванием бит? И почему со знаком минус записывается?

[celen]

Хехе, it's magic.

Надеюсь, мой комментарий ниже про гистограмму убеждает в том, что это *работает*.
Технически, B[offset:offset + size] = -np.log2(rng.random(size)) осуществляет следующую цепочку операций: генерацию случайного числа float64 из равномерного распределения [0;1], взятие от этого логарифма и приведение этого логарифма к целому uint8 (с минусом) через отбрасывание дробной части. Это эквивалентно честной игре Бернулли, но гораздо быстрее. Ещё это эквивалентно тому, что бы взять равномерно распределенный float64 из [0;1], оторвать от его записи экспоненциальную часть и сконвертировать её в uint8. Причем так было бы даже ещё быстрее, но в python я не смог быстро придумать, как эту побитовую операцию грамотно осуществить.

[funca]

Было бы лучше сперва показать наивный алгоритм, отражающий формулировку задачи. А уже потом делать оптимизации. Сейчас статья выглядит как набор глав, между которыми пропущены ключевые связки.

[celen]

Когда я публиковал код из jupyter notebook, я это делал не сколько для чтения непосредственно кода, сколько для его проверки на вашей стороне при наличии такого желания. Для этого коду лучше быть быстрым, то есть немного магическим, с чанками и логарифмическими операциями. Статья достаточно сильно описывает логику происходящего текстом и формулами, что бы разобраться в происходящем вообще без анализа кода.

А вообще эта новая фича хабра sourcecraft достойно справляется с задачей объяснения кода.

[funca]

не сколько для чтения непосредственно кода, сколько для его проверки на вашей стороне

Зависит от того, что понимать под проверкой кода. В представленном варианте не понятно как этот код вообще согласуется с исходной задачей. В статье нет таких формул, которые бы доказывали эквивалентность оптимизированной версии исходному алгоритму. Что тут можно проверить?

Смущает ещё то, что сам по себе питон не ограничивает разрядность целых и позволяет работать с бесконечными последовательностями. В самой задаче вроде нет ни чего, что бы вынуждало держать в памяти результаты всех игр. Поэтому описанние возникших проблем с памятью и разряднотстью в numpy выглядит как какой курьез не по делу. Можно-ли было обойтись вообще без numpy?

[sic]

Ну считайте что я все проверил. У меня по-другому сделано, но результаты совпадают. Вплоть до того, что небольшой недочет и у меня и автора одинаковый, больше 2^64 выигрыши не учитывются. Но, у меня уже больше часа программа работает и только 2^41 выбилось и вот уже триллиарды игр. До 2^64 я все равно не досчитаю никогда (порядка сотен лет с такой скоростью, да и смысла никакого в этом нет, log2(T)/2 там и есть).

А в памяти нужно держать из-за того что автор по этим данным несколько исследований делает, и результаты должны быть согласованны. Кроме того питон не ахти как быстро данные игр генерит даже после его оптимизаций, а наивный подход был бы еще на порядки медленее. И если LLM спросить, про трюк с логарифмом она расскажет.

[celen]

2^41? Мое уважение.

[celen]

Мой оптимизированный код для вычисления B может сбивать с толку, но он дает корректную гистограмму B. Вот какое распределение я получаю :

>> print(B_values)
array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16,
       17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32],
      dtype=uint8)
>> print(B_counts/T)
array([4.99996199e-01, 2.50013790e-01, 1.24990516e-01, 6.25037584e-02,
       3.12467129e-02, 1.56276303e-02, 7.81149906e-03, 3.90555570e-03,
       1.95188215e-03, 9.76342708e-04, 4.87847719e-04, 2.44213734e-04,
       1.22004421e-04, 6.09667040e-05, 3.05534340e-05, 1.52424909e-05,
       7.60983676e-06, 3.83588485e-06, 1.91642903e-06, 9.79984179e-07,
       4.75673005e-07, 2.36090273e-07, 1.16880983e-07, 5.98374754e-08,
       2.56113708e-08, 1.39698386e-08, 8.38190317e-09, 3.72529030e-09,
       1.62981451e-09, 6.98491931e-10, 2.32830644e-10, 2.32830644e-10,
       2.32830644e-10])

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

[strikingamaranth]

Пока толстый сохнет, худой сдохнет -- Старинная русская поговорка

вроде бы все эти поговорки были записаны в конце 19 века на нижегородской ярмарке со слов нетрезвых проституток и в обиходе до этого не были