Como posso fazer uma amostra de uma distribuição com CDF incomputável?

8

Problema relacionado à simulação de semi-informática aqui.

Eu tenho uma distribuição onde

P (x) = $\frac{(e^b-1) e^{b (n-x)}}{e^{b n+b}-1}$

para algumas constantes b e n, e x é um número inteiro tal que $0\leq x \leq n$ .

Agora, preciso provar dessa distribuição. Ele tem um CDF invertível, portanto é possível fazer isso diretamente na teoria. O problema é que os números envolvidos são GRANDES. Tão grande, de fato, que ambas transbordam variáveis formatadas convencionalmente e levam pelo menos minutos (em algum momento eu desisti ...) para calcular usando formatos de precisão arbitrários. Basicamente, o CDF inverso ainda envolve um termo de $e^{b(n+1)}$ , para $350 < n < 3500$ . Apesar disso, os números de saída ainda estarão no intervalo $0-n$ , portanto parece que deve haver uma maneira de fazer isso.

O que estou procurando é uma maneira de amostrar aproximadamente dessa distribuição que é computável. Existem métodos alternativos de amostragem? O que eles são?

sampling simulation monte-carlo John Doucette
fonte

2

Você já pensou em normalizar ou dimensionar seus dados de alguma forma para reduzir o domínio?

EngrStudent

7

O CDF é prontamente invertível. Uma fórmula para a inversão leva ao que deve ser uma das soluções mais simples e convenientes possíveis.

$k$ $0 \le k \le n$ $e^{-b k}$ $q$ $0$ $q_{\max}=\sum_{k=0}^{n} e^{-b k}$ $(1 - e^{-b(n+1)})/(1 - e^{-b})$ $k$

q \geq \sum_{i = 0}^{k} e^{- b i} = \frac{1 - e^{- (k + 1) b}}{1 - e^{- b}} .

$q \ge \sum_{i=0}^{k} e^{-bi} = \frac{1 - e^{-(k+1)b}}{1 -e^{-b}}.$

Álgebra simples fornece a solução

k = - ceiling (\frac{\log (1 - q (1 - e^{- b}))}{b}) .

$k = - \text{ceiling}\left(\frac{\log(1 - q (1-e^{-b}))}{b}\right).$

Aqui está uma Rimplementação construída como todos os outros geradores de números aleatórios: seu primeiro argumento especifica quantos valores de iid a serem gerados e o restante dos argumentos nomeia os parâmetros ( as e as ): $b$ b $n$ n.max

rgeom.truncated <- function(n=1, b, n.max) {
  a <- 1 - exp(-b)
  q.max <- (1 - exp(-b*(n.max+1))) / a
  q <- runif(n, 0, q.max)
  return(-ceiling(log(1 - q*a) / b))
}

Como exemplo de seu uso, vamos gerar um milhão de variáveis aleatórias de acordo com esta distribuição:

b <- 0.001
n.max <- 3500
n.sim <- 10^6
set.seed(17)
system.time(sim <- rgeom.truncated(n.sim, b,n.max))

( Foram necessários segundos.) $0.10$

h <- hist(sim+1, probability=TRUE, breaks=50, xlab="Outcome+1")
pmf <- exp(-b * (0: n.max)); pmf <- pmf / sum(pmf)
lines(0:n.max, pmf, col="Red", lwd=2)

Histograma

( foi adicionado a cada valor para criar um histograma melhor: o procedimento tem uma idiossincrasia (= bug) na qual a primeira barra é muito alta quando o ponto de extremidade esquerdo é definido como zero.) A curva vermelha é a distribuição de referência que esta simulação tenta se reproduzir. Vamos avaliar a qualidade do ajuste com um teste do qui-quadrado: $1$ Rhist

observed <- table(sim)
expected <- n.sim * pmf
chi.square <- (observed-expected)^2 / expected
pchisq(sum(chi.square), n.max, lower.tail=FALSE)

O valor de p é : um ajuste bonito. $0.84$

whuber
fonte

3

Ótima solução. Eu nunca soube que alguém poderia provar dessa maneira (ou seja, baseando-se em amostras de vez de ), mas isso é óbvio em retrospecto.

U n i (0, k), k > 1

$Uni(0, k), k>1$

U n i (0, 1)

$Uni(0,1)$

precisa saber é o seguinte

6

Você está lidando com uma distribuição geométrica truncada com . Existem várias maneiras de abordar isso. $p = 1-e^{-b}$

Eu aconselharia diferentes opções em diferentes situações; algumas opções envolveriam a simulação de um ponto geométrico e a regeneração quando estiver fora do intervalo, levando a parte inteira de um exponencial truncado apropriado ( como aqui ) ou usando qualquer uma das várias técnicas rápidas personalizadas para distribuições discretas em um intervalo finito. Dado que é grande, tomar o piso de um exponencial truncado provavelmente será relativamente rápido, mas se é a melhor opção também depende de . $n$ $b$

Aqui está uma pergunta relacionada sobre math.SE

Antes de tentar sugestões específicas, qual é o intervalo típico de valores para ? $b$

Glen_b -Reinstate Monica
fonte

Obrigado pela sua resposta! b ~ ln (1 + epsilon), onde epsilon é um parâmetro adicional> 0.

John Doucette

1

Então, você converteu minha pergunta sobre o intervalo típico de b em um sobre o intervalo típico de ε. Antes de tentar sugestões específicas, qual é o intervalo típico de valores para ε?

Glen_b -Reinstala Monica

A razão pela qual pergunto é quais abordagens específicas são mais eficientes depende das características da situação. Parece que você já está feliz com a outra resposta, então talvez não valha a pena se preocupar com uma eficiência potencial adicional.

Glen_b -Reinstala Monica

@ JohnDoucette: Se b é quase zero, sua distribuição é quase uniforme em portanto, você pode usar o uniforme como uma proposta em um algoritmo de aceitação de rejeição, pois o limite superior não deve ser terrível.

{0, \dots, n \]

$\{0,\ldots,n\]$

Xi'an

1

@ Xi'an Você precisaria de razoavelmente pequeno em vez de antes que fosse apropriado usar uma distribuição uniforme, porque a taxa de aceitação é , que será ineficientemente baixo quando .

n b

$nb$

b \approx 0

$b\approx 0$

(1 - e^{- (n + 1) b}) / ((n + 1) (1 - e^{- b}))

$(1- e^{-(n+1)b})/((n+1)(1-e^{-b}))$

\approx (1 - \exp (- n b)) / (n b)

$\approx (1-\exp(-n b))/(n b)$

n b ≫ 1

$nb \gg 1$

whuber

4

Primeiro, observe que que, se fosse contínuo, estaria relacionado a uma distribuição exponencial. Então, o que você pode fazer é simular a partir de uma distribuição exponencial truncada e pegar a (parte inteira) das observações. $P(x)\propto e^{-bx}$ $x$ floor()

O cdf de um exponencial truncado é

F (x; n, b) = \frac{1 - e^{- b x}}{1 - e^{- b n}} .

$F(x;n,b)= \dfrac{1-e^{-bx}}{1-e^{-bn}}.$

Então, se fizermos , obtemos que . Se for grande, que sugere aproximar . $F(x;n,b)=u$ $x=-\dfrac{1}{b}\log[1-u(1-e^{-bn})]$ $bn$ $e^{-bn}\approx 0$ $x\approx -\dfrac{1}{b}\log[1-u]$

rweirdp <- function(ns,n,b){
u <- runif(ns)
samp <- - log(1-u*(1-exp(-n*b)))/b
return(floor(samp))
}

rweirdp(1000,10,1)

Pessoa
fonte

Eu acho que isso é basicamente o que eu estava procurando. bn é sempre muito grande, amostragem proporcional faz sentido. Não estava ciente do mapeamento, embora seja claro em retrospecto. Obrigado!

John John Doucette

Fico feliz em ver que isso ajudou. Acho que não expliquei corretamente, mas essa abordagem produz amostras a partir da distribuição exata do alvo. Felicidades.

Pessoa

@ Xi'an Os pesos não são os mesmos se alguém usar o valor e pegar a parte inteira?

e^{- b n}

$e^{-bn}$

Pessoa

@ Xi'an Acho que aparece prazo no numerador de , até um factorisation ...

P (x)

$P(x)$

Person

1

@ Xi'an Na verdade, este trabalho fornecido rweirdpé modificado para mudar npara n+1. (Como indicado aqui, ele nunca retornará um valor igual a n: esse é o efeito da aproximação.) Uma análise um pouco mais rigorosa é dada em minha resposta. Embora eu obtenha uma fórmula com aparência diferente, ela é equivalente à (mais simples!), Dada aqui, uma vez que a modificação n-> n+1é feita.

whuber

4

Uma maneira de obter amostra da distribuição de destino é $p(k)\propto \exp\{-bk\}$

execute um experimento Metropolis-Hastings para determinar o suporte (interessante) da distribuição, ou seja, em qual subconjunto de ele se concentra; $\{0,1,\ldots,n\}$

metro=function(N,b,n){
x=sample(0:n,N,rep=TRUE)
for (t in 2:N){
  x[t]=prop=x[t-1]+sample(c(-1,1),1)

  if ((prop<0)||(prop>n)||(log(runif(1))>b*(x[t]-prop)))
      x[t]=x[t-1]
  }
return(x)
}

Use o suporte assim determinado, digamos, para calcular as probabilidades exatas como para evitar estouros. $\{k_0,\ldots,k_1\}$ $p(k)\propto \exp\{-bk+bk_0\}$

Atualização: Ao pensar mais sobre isso, como está diminuindo em k, o suporte efetivo da distribuição sempre começará em . Se for muito grande, esse suporte terminará muito rapidamente, caso em que não importa muito, pois grandes valores de nunca serão visitados. Se for muito pequeno, o pdf será quase plano, o que significa que é possível usar uma distribuição uniforme em como uma proposta de aceitação / rejeição. E use logs na etapa de aceitação para evitar estouros. $p(\cdot)$ $k_0=0$ $b$ $n$ $k$ $b$ $\{0,1,\ldots,n\}$

Xi'an
fonte

Como posso fazer uma amostra de uma distribuição com CDF incomputável?

Respostas: