Simule uma variável de Bernoulli com probabilidade usando uma moeda polarizada

9

Alguém pode me dizer como simular , onde , usando um sorteio (quantas vezes você precisar) com ?Bernoulli(ab)a,bNP(H)=p

Eu estava pensando em usar amostras de rejeição, mas não consegui identificá-las.

Abracadabra
fonte
11
Esta é uma pergunta originalmente de um curso ou livro? Em caso afirmativo, adicione a [self-study]tag e leia seu wiki . Observe que não há necessidade de pedir ajuda no final de sua pergunta - sabemos que todos os que postam aqui estão esperando por ajuda!
Silverfish
11
Há um excelente post de @Glen_b aqui em algum lugar (embora eu não me lembro onde) sobre por que não existe uma "moeda tendenciosa com probabilidade ", mas eu sei que isso é apenas uma questão periférica para sua pergunta! p
Silverfish
2
@dsaxton A pergunta diz "quantos você precisar"; será finito com probabilidade 1, mas não limitado (você pode exceder qualquer número fixo de lançamentos), mas objetar com base nisso seria como dizer "jogar uma moeda justa até que você receba uma cabeça" não é viável como um método de geração geométrica ( números aleatórios.12
Glen_b -Reinstata Monica 25/04
11
@AbracaDabra Este é um exercício para uma aula? Se não, como surge?
Glen_b -Reinstala Monica
11
@Glen_b: Não é um exercício da minha turma. Apenas me ocorreu nesta cadeia de pensamento ...: De acordo com a probabilidade clássica, pegue uma moeda justa, à medida que aumenta o número de jogadas, a proporção de converge para a metade. Portanto, também deve ser verdade para os tendenciosos ... Isso significa que uma moeda converge para um número específico, você precisa que seja esse número. Agora pensei: e se quisermos produzir um número, mas tivermos uma moeda com outro número (conhecido ou desconhecido)? #Heads#tailsP(H)P(H)
AbracaDabra

Respostas:

8

Como existem inúmeras soluções, vamos encontrar uma eficiente .

A idéia por trás disso começa com uma maneira padrão de implementar uma variável de Bernoulli: compare uma variável aleatória uniforme com o parâmetro . Quando , retorne ; caso contrário, retorne .Ua/bU<a/b10

Podemos usar a moeda como um gerador uniforme de números aleatóriosp . Para gerar um número uniformemente dentro de qualquer intervalo , jogue a moeda. Quando estiver no cabeçalho, gere recursivamente um valor uniforme na primeira parte do intervalo; quando for coroa, gere recursivamente a partir da última parte de do intervalo. Em algum momento, o intervalo de destino se tornará tão pequeno que realmente não importa como você escolhe um número: é assim que a recursão começa. É óbvio que esse procedimento gera variáveis ​​uniformes (até a precisão desejada), como é facilmente comprovado por indução.U[x,y)XpX1p

Essa ideia não é eficiente, mas leva a um método eficiente. Como em cada estágio você desenha um número de um determinado intervalo , por que não verificar primeiro se é necessário desenhá-lo? Se o seu valor-alvo estiver fora desse intervalo, você já saberá o resultado da comparação entre o valor aleatório e o alvo. Assim, esse algoritmo tende a terminar rapidamente. (Isso pode ser interpretado como o procedimento de amostragem por rejeição solicitado na pergunta.)[x,y)

Podemos otimizar ainda mais esse algoritmo. Em qualquer estágio, na verdade, temos duas moedas que podemos usar: ao re-rotular nossa moeda, podemos transformá-la em uma que tem cara de chance . Portanto, como precomputação, podemos escolher recursivamente qualquer remarcação que leve ao menor número esperado de inversões necessárias para a rescisão. (Este cálculo pode ser uma etapa cara.)1p

Por exemplo, é ineficiente usar uma moeda com para emular diretamente uma variável de Bernoulli : são necessários quase dez movimentos em média. Mas se usarmos uma moeda , em apenas dois lançamentos teremos certeza de que será feito e o número esperado de lançamentos será de apenas .p=0.9(0.01)p=10.0=0.11.2

Aqui estão os detalhes.

Particionar qualquer intervalo semi-aberto nos intervalosI=[x,y)

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

Isto define as duas transformações es(,H)s(,T) que operam em intervalos de semi-abertas.

I

t<I

tIt<xxIt>ItI

a/b=tt0t1

Z

  1. n=0In=I0=[0,1)

  2. Enquanto {Jogue a moeda para produzir . Defina Incremento .}(tIn)Xn+1In+1=S(In,Xn+1).n

  3. Se , defina . Caso contrário, defina .t>In+1Z=1Z=0


Implementação

Para ilustrar, aqui está uma Rimplementação do aloritmo como a função draw. Seus argumentos são o valor alvo o intervalo , inicialmente . Ele usa a função auxiliar implementando . Embora não seja necessário, também rastreia o número de lançamentos de moedas. Retorna a variável aleatória, a contagem de lançamentos e o último intervalo inspecionado.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))
}

Como um exemplo da sua utilização e ensaio da sua exactidão, assumir o caso e . Vamos desenhar valores usando o algoritmo, relatar a média (e seu erro padrão) e indicar o número médio de inversões usadas.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.00950.010.51549.8861p0.00941.177

whuber
fonte
Não consigo deixar de ver as semelhanças entre esta solução e a Solução 2 na minha resposta. Enquanto eu assumo uma moeda imparcial (PS, solução realmente interessante para o problema da moeda tendenciosa) e faço todos os cálculos / comparações na base 2, você faz todos os cálculos / comparações na base 10. Quais são seus pensamentos?
precisa saber é o seguinte
11
@cam Eu acho que você pode ser enganado pelos meus exemplos: embora eles usem bons números na base 10, a construção não tem nada a ver com uma base específica.
whuber
2
a/bpn(1p)m(n+mm)pn(1p)m
4

Aqui está uma solução (um pouco bagunçada, mas é minha primeira facada). Você pode realmente ignorar WLOG assumem . Por quê? Existe um algoritmo inteligente para gerar um lançamento de moeda imparcial a partir de dois lançamentos de moedas tendenciosos. Portanto, podemos assumir .P(H)=pP(H)=1/2P(H)=1/2

Para gerar um , posso pensar em duas soluções (a primeira não é minha, mas a segunda é uma generalização):Bernoulli(ab)

Solução 1

Vire a moeda imparcial vezes. Se cabeça não estiver presente, comece de novo. Se cabeças estão presentes, o retorno se a primeira moeda é uma ou cabeças não (porque )baaP(first coin is heads | a heads in b coins)=ab

Solução 2

Isso pode ser estendido para qualquer valor de . Escreva em formato binário. Por exemplo,Bernoulli(p)p0.1=0.0001100110011001100110011...base 2

Criaremos um novo número binário usando lançamentos de moedas. Comece com e adicione dígitos, dependendo se uma cara (1) ou coroa (0) aparecer. Em cada flip, compare seu novo número binário com a representação binária de até o mesmo dígito . Eventualmente, os dois divergem e retornam se o for maior que o seu número binário.0.p bin(p)

Em 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

Alguma prova:

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

é de cerca de 0,4 (no entanto, não é rápido)

Cam.Davidson.Pilon
fonte
Boa resposta, mas você pode explicar com seu método 1 como fazer por p irracional?
AbracaDabra
2
@AbracaDabra, por que a racionalidade de importa? p
Glen_b -Reinstala Monica
@ AbracaDabra: qualquer que seja o valor de , a probabilidade de obter e é a mesma, ou seja, , portanto, a probabilidade de obter um contra o outro é . p(0,1)(1,0)p(1p)1/2
Xian
4

Vejo uma solução simples, mas sem dúvida há muitas maneiras de fazê-lo, algumas presumivelmente mais simples que isso. Essa abordagem pode ser dividida em duas etapas:

  1. Gerando a partir de dois eventos com igual probabilidade, dado um procedimento injusto de lançamento de moeda (a combinação da moeda específica e o método pelo qual ela é lançada, gerando uma cabeça com probabilidade ). Podemos chamar esses dois eventos igualmente prováveis e . [Existe uma abordagem simples para isso, que exige fazer pares de arremessos e para produzir dois resultados igualmente prováveis, com todos os outros resultados levando à geração de um novo par de rolos para tentar novamente.]pHTH=(H,T)T=(T,H)

  2. Agora você gera uma caminhada aleatória com dois estados absorventes usando a moeda justa simulada. Ao escolher a distância dos estados de absorção da origem (um acima e outro abaixo), você pode definir a chance de absorção, digamos que o estado de absorção superior seja a proporção desejada de números inteiros. Especificamente, se você colocar a barreira de absorção superior em e a inferior em (e iniciar o processo a partir da origem) e executar a caminhada aleatória até a absorção, a probabilidade de absorção na barreira superior é .a(ba)aa+(ba)=ab

    (Existem alguns cálculos a serem feitos aqui para mostrá-lo, mas você pode obter facilmente as probabilidades trabalhando com relações de recorrência ... ou pode fazê-lo somando séries infinitas ... ou de outras maneiras.)

Glen_b -Reinstate Monica
fonte