Хабрахабр

[Перевод] Как выбрать случайное число от 1 до 10

Представьте, что вам нужно сгенерировать равномерно распределённое случайное число от 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, также именуется линейной оптимизацией) — метод достижения наилучшего результата… в математической модели, требования которой представлены линейными отношениями… Стандартная форма представляет собой обычную и наиболее интуитивную форму описания задачи линейного программирования.

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

Аналогично можно сформулировать и проблему перераспределения.
У нас есть набор переменных $(x_$, каждая из которых кодирует долю вероятности, перераспределённую от целого числа $i$ (от 1 до 10) к целому числу $j$ (от 1 до 10). Поэтому, если $(x_{7,1} = 0.2$, то нам нужно перенести 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%. Другими словами, для каждого $j$ от 1 до 10:

1$" data-tex="display"/> <img src="https://habrastorage.org/getpro/habr/formulas/266/e9e/845/266e9e84523992398c85afdb20ef9ece.svg" alt="$x_{1, j} + x_{2, j} + \ldots\ x_{10, j} = 0.

Можем представить эти ограничения в виде списка массивов в 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) })

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

1$" data-tex="display"/> <img src="https://habrastorage.org/getpro/habr/formulas/266/e9e/845/266e9e84523992398c85afdb20ef9ece.svg" alt="$x_{1, j} + x_{2, j} + \ldots\ x_{10, j} = 0.

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 (x_{1, 1} + x_{2, 2} + \ldots\ x_{10, 10})$

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")

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

Но у нас только комната, полная людей. Было бы довольно просто решить задачу, если бы у нас был доступ к генератору равномерно распределённых случайных чисел (от 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!

Теги
Показать больше

Похожие статьи

Добавить комментарий

Ваш e-mail не будет опубликован. Обязательные поля помечены *

Кнопка «Наверх»
Закрыть