Как найти произвольное число

Время на прочтение
6 мин

Количество просмотров 35K

Представьте, что вам нужно сгенерировать равномерно распределённое случайное число от 1 до 10. То есть целое число от 1 до 10 включительно, с равной вероятностью (10%) появления каждого. Но, скажем, без доступа к монетам, компьютерам, радиоактивному материалу или другим подобным источникам (псевдо) случайных чисел. У вас есть только комната с людьми.

Предположим, что в этой комнате чуть более 8500 студентов.

Самое простое — попросить кого-нибудь: «Эй, выбери случайное число от одного до десяти!». Человек отвечает: «Семь!». Отлично! Теперь у вас есть число. Однако вы начинаете задаваться вопросом, является ли оно равномерно распределённым?

Поэтому вы решили спросить ещё несколько человек. Вы продолжаете их спрашивать и считать их ответы, округляя дробные числа и игнорируя тех, кто думает, что диапазон от 1 до 10 включает 0. В конце концов вы начинаете видеть, что распределение вообще не равномерное:

library(tidyverse)

probabilities <-
  read_csv("https://git.io/fjoZ2") %>%
  count(outcome = round(pick_a_random_number_from_1_10)) %>%
  filter(!is.na(outcome),
         outcome != 0) %>%
  mutate(p = n / sum(n))

probabilities %>%
  ggplot(aes(x = outcome, y = p)) +
  geom_col(aes(fill = as.factor(outcome))) +
  scale_x_continuous(breaks = 1:10) +
  scale_y_continuous(labels = scales::percent_format(),
                     breaks = seq(0, 1, 0.05)) +
  scale_fill_discrete(h = c(120, 360)) +
  theme_minimal(base_family = "Roboto") +
  theme(legend.position = "none",
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank()) +
  labs(title = '"Pick a random number from 1-10"',
       subtitle = "Human RNG distribution",
       x = "",
       y = NULL,
       caption = "Source: https://www.reddit.com/r/dataisbeautiful/comments/acow6y/asking_over_8500_students_to_pick_a_random_number/")

Данные с Reddit

Вы хлопаете себя по лбу. Ну конечно, оно не будет случайным. В конце концов, нельзя доверять людям.

Итак, что делать?

Вот бы найти какую-то функцию, которая преобразует распределение «человеческого ГСЧ» в равномерное распределение…

Интуиция тут относительно проста. Нужно всего лишь взять массу распределения оттуда, где она выше 10%, и переместить туда, где она меньше 10%. Так, чтобы все столбцы на графике были одного уровня:

По идее, такая функция должна существовать. Фактически, должно быть много различных функций (для перестановки). В крайнем случае, можно «разрезать» каждый столбец на бесконечно малые блоки и построить распределение любой формы (как кирпичики Lego).

Конечно, такой экстремальный пример немного громоздок. В идеале мы хотим сохранить как можно больше исходного распределения (т. е. сделать как можно меньше измельчений и перемещений).

Как найти такую функцию?

Ну, наше объяснение выше звучит очень похоже на линейное программирование. Из Википедии:

Линейное программирование (LP, также именуется линейной оптимизацией) — метод достижения наилучшего результата… в математической модели, требования которой представлены линейными отношениями… Стандартная форма представляет собой обычную и наиболее интуитивную форму описания задачи линейного программирования. Она состоит из трёх частей:

• Линейная функция, которую необходимо максимизировать
• Проблемные ограничения следующей формы
• Неотрицательные переменные

Аналогично можно сформулировать и проблему перераспределения.

Представление проблемы

У нас есть набор переменных

, каждая из которых кодирует долю вероятности, перераспределённую от целого числа

(от 1 до 10) к целому числу

(от 1 до 10). Поэтому, если

, то нам нужно перенести 20% ответов от семёрки к единице.

variables <-
  crossing(from = probabilities$outcome,
           to   = probabilities$outcome) %>%
  mutate(name = glue::glue("x({from},{to})"),
         ix = row_number())

variables

## # A tibble: 100 x 4
##     from    to name       ix
##    <dbl> <dbl> <glue>  <int>
##  1     1     1 x(1,1)      1
##  2     1     2 x(1,2)      2
##  3     1     3 x(1,3)      3
##  4     1     4 x(1,4)      4
##  5     1     5 x(1,5)      5
##  6     1     6 x(1,6)      6
##  7     1     7 x(1,7)      7
##  8     1     8 x(1,8)      8
##  9     1     9 x(1,9)      9
## 10     1    10 x(1,10)    10
## # … with 90 more rows

Мы хотим ограничить эти переменные таким образом, чтобы все перераспределённые вероятности суммировались в 10%. Другими словами, для каждого

от 1 до 10:

Можем представить эти ограничения в виде списка массивов в R. Позже свяжем их в матрицу.

fill_array <- function(indices,
                       weights,
                       dimensions = c(1, max(variables$ix))) {
  init <- array(0, dim = dimensions)

  if (length(weights) == 1) {
    weights <- rep_len(1, length(indices))
  }

  reduce2(indices, weights, function(a, i, v) {
    a[1, i] <- v
    a
  }, .init = init)
}

constrain_uniform_output <-
  probabilities %>%
  pmap(function(outcome, p, ...) {
    x <-
      variables %>%
      filter(to == outcome) %>%
      left_join(probabilities, by = c("from" = "outcome"))

    fill_array(x$ix, x$p)
  })

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

в диапазоне от 1 до 10:

one_hot <- partial(fill_array, weights = 1)

constrain_original_conserved <-
  probabilities %>%
  pmap(function(outcome, p, ...) {
    variables %>%
      filter(from == outcome) %>%
      pull(ix) %>%
      one_hot()
  })

Как уже говорилось, мы хотим максимизировать сохранение исходного распределения. Это наша цель (objective):

maximise_original_distribution_reuse <-
  probabilities %>%
  pmap(function(outcome, p, ...) {
    variables %>%
      filter(from == outcome,
             to == outcome) %>%
      pull(ix) %>%
      one_hot()
  })

objective <- do.call(rbind, maximise_original_distribution_reuse) %>% colSums()

Затем передаём проблему солверу, например, пакету lpSolve в R, объединив созданные ограничения в одну матрицу:

# Make results reproducible...
set.seed(23756434)

solved <- lpSolve::lp(
  direction    = "max",
  objective.in = objective,
  const.mat    = do.call(rbind, c(constrain_original_conserved, constrain_uniform_output)),
  const.dir    = c(rep_len("==", length(constrain_original_conserved)),
                   rep_len("==", length(constrain_uniform_output))),
  const.rhs    = c(rep_len(1, length(constrain_original_conserved)),
                   rep_len(1 / nrow(probabilities), length(constrain_uniform_output)))
)

balanced_probabilities <-
  variables %>%
  mutate(p = solved$solution) %>%
  left_join(probabilities,
            by = c("from" = "outcome"),
            suffix = c("_redistributed", "_original"))

Возвращается следующее перераспределение:

library(gganimate)

redistribute_anim <-
  bind_rows(balanced_probabilities %>%
            mutate(key   = from,
                   state = "Before"),
            balanced_probabilities %>%
            mutate(key   = to,
                   state = "After")) %>%
  ggplot(aes(x = key, y = p_redistributed * p_original)) +
  geom_col(aes(fill = as.factor(from)),
           position = position_stack()) +
  scale_x_continuous(breaks = 1:10) +
  scale_y_continuous(labels = scales::percent_format(),
                     breaks = seq(0, 1, 0.05)) +
  scale_fill_discrete(h = c(120, 360)) +
  theme_minimal(base_family = "Roboto") +
  theme(legend.position = "none",
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank()) +
  labs(title = 'Balancing the "Human RNG distribution"',
       subtitle = "{closest_state}",
       x = "",
       y = NULL) +
  transition_states(
    state,
    transition_length = 4,
    state_length = 3
  ) +
  ease_aes('cubic-in-out')

animate(
  redistribute_anim,
  start_pause = 8,
  end_pause = 8
)

Отлично! Теперь у нас есть функция перераспределения. Давайте поближе посмотрим, как именно движется масса:

balanced_probabilities %>%
  ggplot(aes(x = from, y = to)) +
  geom_tile(aes(alpha = p_redistributed, fill = as.factor(from))) +
  geom_text(aes(label = ifelse(p_redistributed == 0, "", scales::percent(p_redistributed, 2)))) +
  scale_alpha_continuous(limits = c(0, 1), range = c(0, 1)) +
  scale_fill_discrete(h = c(120, 360)) +
  scale_x_continuous(breaks = 1:10) +
  scale_y_continuous(breaks = 1:10) +
  theme_minimal(base_family = "Roboto") +
  theme(panel.grid.minor = element_blank(),
        panel.grid.major = element_line(linetype = "dotted"),
        legend.position = "none") +
  labs(title = "Probability mass redistribution",
       x = "Original number",
       y = "Redistributed number")

Эта диаграмма говорит, что примерно в 8% случаев, когда кто-то называет восемь в качестве случайного числа, вам нужно воспринимать ответ как единицу. В остальных 92% случаев он остаётся восьмёркой.

Было бы довольно просто решить задачу, если бы у нас был доступ к генератору равномерно распределённых случайных чисел (от 0 до 1). Но у нас только комната, полная людей. К счастью, если вы готовы примириться с несколькими небольшими неточностями, то из людей можно сделать довольно хороший ГСЧ, не спрашивая более двух раз.

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

probabilities %>%
  transmute(number = outcome,
            probability = scales::percent(p))

## # A tibble: 10 x 2
##    number probability
##     <dbl> <chr>
##  1      1 3.4%
##  2      2 8.5%
##  3      3 10.0%
##  4      4 9.7%
##  5      5 12.2%
##  6      6 9.8%
##  7      7 28.1%
##  8      8 10.9%
##  9      9 5.4%
## 10     10 1.9%

Например, когда кто-то даёт нам восемь в качестве случайного числа, нужно определить, должна ли эта восьмёрка стать единицей или нет (вероятность 8%). Если мы спросим другого человека о случайном числе, то с вероятностью 8,5% он ответит «два». Так что если это второе число равно 2, мы знаем, что должны вернуть 1 как равномерно распределённое случайное число.

Распространив эту логику на все числа, получаем следующий алгоритм:

По этому алгоритму вы можете использовать группу людей для получения чего-то близкого к генератору равномерно распределённых случайных чисел от 1 до 10!

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

Но есть одна проблема, любой алгоритм можно предугадать и он начнет зацикливаться, так в 1983г на шоу «Поймай свою удачу», водитель грузовика Майкл Ларсон он выиграл просто блять сумасшедшие деньги просчитав алгоритм случайных числе этой передачи. После этого случая появился сайт random.org, который генерирует по настоящему рандомное число.

Давайте разберем, как он их получает. Получает этот сайт благодаря матушки природы, да-да наша природа которую мы порой блять недооцениваем, потому-что только матушка природа способна получить хаос который мы не можем предугадать, точнее можем, но это будет пиздец как не скоро, т.к. надо будет рассчитать движение каждого атома во вселенной. Сайт random.org генерирует свои основываясь на данных атмосферных шумов.

Окей, нахуя я тут распинался и рассказывал про этот сайт рандом орг?

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

Спасибо за то, что вы с нами.
С любовью, Рителлинг favorite