Имитировать переменную Бернулли с вероятностью

9

Может кто-нибудь сказать мне, как смоделировать , где , используя подбрасывание монеты (столько раз, сколько вам нужно) с ?Bernoulli(ab)a,bNP(H)=p

Я думал об использовании выборки отклонения, но не мог прибить это.

абракадабра
источник
1
Это вопрос из курса или учебника? Если это так, пожалуйста, добавьте [self-study]тег и прочитайте его вики . Обратите внимание, что нет необходимости ставить просьбу о помощи в конце вашего вопроса - мы знаем, что все, кто публикует здесь, надеются на помощь!
Серебряная рыба
1
Где-то здесь есть отличный пост @Glen_b (хотя я не могу вспомнить где) о том, почему не существует такой вещи, как «предвзятая монета с вероятностью p », но я знаю, что это только второстепенная проблема в вашем вопросе!
Серебряная рыба
2
@dsaxton Вопрос говорит: «столько, сколько вам нужно»; оно будет конечным с вероятностью 1, но не ограниченным (вы можете превысить любое фиксированное количество бросков), но возражение на этой основе будет все равно, что сказать «бросить честную монету, пока вы не получите голову», не является жизнеспособным в качестве метода для генерации геометрического ( 1 случайных числа. 12
Glen_b
1
@AbracaDabra Это упражнение для класса? Если нет, то как это возникает?
Glen_b
1
@Glen_b: это не упражнение из моего класса. Это пришло мне в голову в этой цепочке мыслей ...: Согласно классической вероятности, возьмите справедливую монету, поскольку вы увеличиваете количество бросков, отношение сходится к половине. Так что это должно быть верно и для предвзятых ... Это означает, что монета должна сходиться к определенному числу, вам нужно, чтобыP(H)было этим числом. Теперь я подумал, а что если мы хотим произвести число, но у нас есть монета сP(H)другим числом (известным или неизвестным)? #Heads#tailsP(H)P(H)
AbracaDabra

Ответы:

8

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

Идея, лежащая в основе этого, начинается со стандартного способа реализации переменной Бернулли: сравнить равномерную случайную переменную с параметром a / b . Когда U < a / b , вернуть 1 ; в противном случае верните 0 .Ua/bU<a/b10

Мы можем использовать монету как генератор случайных чиселp . Чтобы сгенерировать число равномерно в любом интервале [ x , y ) , подбросьте монету. Когда это головы, рекурсивно генерировать равномерное значение X в первой части p интервала; когда это хвосты, рекурсивно генерировать X из последних 1 - рU[x,y)XpX1pчасть интервала. В какой-то момент целевой интервал станет настолько маленьким, что на самом деле не имеет значения, как вы выбираете из него число: именно так начинается рекурсия. Очевидно, что эта процедура генерирует однородные переменные (с любой желаемой точностью), что легко подтверждается индукцией.

Эта идея неэффективна, но приводит к эффективному методу. Поскольку на каждом этапе вы собираетесь рисовать число из некоторого заданного интервала , почему бы сначала не проверить, нужно ли вообще его рисовать? Если целевое значение лежит за пределами этого интервала, вы уже знаете результат сравнения между случайным значением и целью. Таким образом, этот алгоритм имеет тенденцию быстро завершаться. (Это может быть истолковано как процедура отбора проб , указанная в вопросе.)[x,y)

Мы можем оптимизировать этот алгоритм дальше. На любом этапе у нас на самом деле есть две монеты, которые мы можем использовать: перемаркировав нашу монету, мы можем превратить ее в монету с вероятностью . Следовательно, в качестве предварительного вычисления мы можем рекурсивно выбрать тот, который приведет к перемаркировке, что приведет к меньшему ожидаемому числу бросков, необходимых для завершения. (Этот расчет может быть дорогим шагом.)1p

Например, неэффективно использовать монету с для прямой эмуляции переменной Бернулли ( 0,01 ) : в среднем это занимает почти десять бросков. Но если мы используем монету с p = 1 - 0.0 = 0.1 , то всего за два броска мы обязательно это сделаем, и ожидаемое количество бросков будет всего 1.2 .p=0.9(0.01)p=10.0=0.11.2

Вот подробности.

Разбиение любого заданного полуоткрытого интервала на интервалыI=[x,y)

[x,y)=[x,x+(yx)p)[x+(yx)p,y)=s(I,H)s(I,T).

Это определяет два преобразования и s ( , T ), которые действуют на полуоткрытых интервалах.s(,H)s(,T)

С точки зрения терминологии, если - любое множество действительных чисел, пусть выражениеI

t<I

означает , что является нижней границей для I : т < х для всех х I . Аналогичным образом , т > I означает , что т является верхней границей для ввода .tIt<xxIt>ItI

Напишите . (Фактически, это не будет иметь значения, если t будет действительным, а не рациональным; нам нужно только, чтобы 0 t 1. )a/b=tt0t1

Вот алгоритм для получения переменной с желаемым параметром Бернулли:Z

  1. Установите и I n = I 0 = [ 0 , 1 ) .n=0In=I0=[0,1)

  2. While {Бросить монету, чтобы получить X n + 1 . Установите I n + 1 = S ( I n , X n + 1 ) . Инкремент n .}(tIn)Xn+1In+1=S(In,Xn+1).n

  3. Если тогда установить Z = 1 . В противном случае установите Z = 0 .t>In+1Z=1Z=0


Реализация

Чтобы проиллюстрировать это, здесь приведена Rреализация алгоритма как функции draw. Его аргументами являются целевое значение и интервал [ x , y ) , изначально [ 0 , 1 ) . Он использует вспомогательную функцию, реализующую s . Хотя это и не нужно, но также отслеживает количество брошенных монет. Он возвращает случайную величину, количество бросков и последний проверенный интервал.t[x,y)[0,1)ss

s <- function(x, ab, p) {
  d <- diff(ab) * p
  if (x == 1) c(ab[1], ab[1] + d) else c(ab[1] + d, ab[2])
}
draw <- function(target, p) {
  between <- function(z, ab) prod(z - ab) <= 0
  ab <- c(0,1)
  n <- 0
  while(between(target, ab)) {
    n <- n+1; ab <- s(runif(1) < p, ab, p)
  }
  return(c(target > ab[2], n, ab))
}

В качестве примера его использования и испытания его точности, возьмите случай и р = 0,9 . Обратим 10 , 000 значения , используя алгоритм, отчет о среднем (и ее стандартной ошибки), и указать среднее число перестроек используется.t=1/100p=0.910,000

target <- 0.01
p <- 0.9
set.seed(17)
sim <- replicate(1e4, draw(target, p))

(m <- mean(sim[1, ]))                           # The mean
(m - target) / (sd(sim[1, ]) / sqrt(ncol(sim))) # A Z-score to compare to `target`
mean(sim[2, ])                                  # Average number of flips

В этом моделировании сальто были головы. Несмотря на то, что показатель ниже 0,01 , Z-показатель - 0,5154 не является значимым: это отклонение можно объяснить случайностью. Среднее количество сальто составило 9,886 - чуть меньше десяти. Если бы мы использовали 1 - р монетку, средний было бы 0,0094 --still не значительно отличается от цели, но только 1,177 щелчки были бы необходимы в среднем.0.00950.010.51549.8861p0.00941.177

Whuber
источник
Я не могу не увидеть сходство между этим решением и решением 2 в моем ответе. В то время как я предполагаю несмещенную монету (PS действительно интересное решение проблемы смещенной монеты), и делаю все вычисления / сравнения в базе-2, вы делаете все вычисления / сравнения в базе 10. Что вы думаете?
Cam.Davidson.Pilon
1
@cam Я думаю, что вы можете быть обмануты моими примерами: хотя они используют хорошие числа в базе 10, конструкция не имеет никакого отношения к какой-либо конкретной базе.
whuber
2
(+1) Очень аккуратное разрешение. Оптимизация заключается в верхнем и нижнем ограничении такими степенями, как p n ( 1 - p ) m и / или ( n + m).a/bpn(1p)m. Было бы неплохо найти оптимальную дихотомию с точки зрения количества смоделированных Бернулли. (n+mm)pn(1p)m
Сиань
5

Вот решение (немного грязное, но это мой первый удар). Вы можете фактически игнорировать и без потери общности , предположим , Р ( Н ) = 1 / 2 . Почему? Существует умный алгоритм для генерации несмещенного подбрасывания монет из двух предвзятых подбрасываний монет. Таким образом , мы можем предположить , P ( H ) = 1 / 2 .P(H)=pP(H)=1/2P(H)=1/2

Для того, чтобы сгенерировать я могу придумать два решения (первое не мое, а второе обобщение):Bernoulli(ab)

Решение 1

Переверните несмещенную монету раз. Если а голова нет, начните заново. Если через головку может присутствовать, возвратом первой монетой является ли головами или нет (потому что Р ( первая монета головы |  а  головы в  б  монетах ) =baa )P(first coin is heads | a heads in b coins)=ab

Решение 2

Это может быть распространено на любое значение . Напишите p в двоичной форме. Например, 0,1 = 0,0001100110011001100110011 ... основание 2Bernoulli(p)p0.1=0.0001100110011001100110011...base 2

Мы создадим новое двоичное число, используя броски монет. Начните с и добавьте цифры в зависимости от того, появляются ли головы (1) или хвосты (0). При каждом броске сравнивайте новое двоичное число с двоичным представлением p до той же цифры . В конце концов они будут расходиться и возвращаться, если b i n ( p ) больше вашего двоичного числа.0.p bin(p)

В Python:

def simulate(p):
    binary_p = float_to_binary(p)
    binary_string = '0.'
    index = 3
    while True:
        binary_string += '0' if random.random() < 0.5 else '1'
        if binary_string != binary_p[:index]:
            return binary_string < binary_p[:index]
        index += 1

Некоторые доказательства:

np.mean([simulate(0.4) for i in range(10000)])

около 0,4 (не быстро, однако)

Cam.Davidson.Pilon
источник
Хороший ответ, но вы можете объяснить своим методом 1, как это сделать для иррационального p?
AbracaDabra
2
@AbracaDabra, почему рациональность имеет значение? p
Glen_b
@AbracaDabra: независимо от значения , вероятность получения ( 0 , 1 ) и ( 1 , 0 ) являются такими же, а именно : р ( 1 - р ) , следовательно , вероятность получения друг против друга на 1 / 2 . p(0,1)(1,0)p(1p)1/2
Сиань
4

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

  1. pHTH=(H,T)T=(T,H)

  2. a(ba)aa+(ba)=ab

    (Здесь нужно сделать некоторые вычисления, чтобы показать это, но вы можете довольно легко вывести вероятности, работая с рекуррентными соотношениями ... или вы можете сделать это, суммируя бесконечные ряды ... или есть другие способы.)

Glen_b - Восстановить Монику
источник