Como amostrar de

Quero amostrar de acordo com uma densidade

Alguém sabe como tirar amostras dessa densidade facilmente? Talvez seja padrão e apenas algo que eu não conheço?

Eu posso pensar em um algoritmo de amostragem de rejeição estúpido que funcionará mais ou menos (encontre o modo $a^*$ de $f$ , amostra $(a,u)$ de uniforme em uma caixa grande $[0,10a^*]\times [0,f(a^*)]$ e rejeite se $u>f(a)$ ), mas (i) não é de todo eficiente e (ii) $f(a^*)$será muito grande para um computador suportar facilmente até mesmo moderadamente grandes $c$ e $d$ . (Note-se que o modo para a grande $c$ e $d$ é, aproximadamente, a $a=cd$ ).

Agradecemos antecipadamente por qualquer ajuda!

A amostragem por rejeição funcionará excepcionalmente bem quando e é razoável para c d ≥ $c d \ge \exp(5)$ .$c d \ge \exp(2)$

Para simplificar um pouco a matemática, deixe , escreva x = a e observe que$k = c d$$x = a$

para . Definindo x = L 3 / 2 dá$x \ge 1$$x = u^{3/2}$

para . Quando k ≥ exp ( 5 ) , essa distribuição é extremamente próxima de Normal (e se aproxima à medida que k aumenta). Especificamente, você pode$u \ge 1$$k \ge \exp(5)$$k$

1. Encontre o modo de numericamente (usando, por exemplo, Newton-Raphson).$f(u)$

2. Expanda para a segunda ordem sobre seu modo.$\log{f(u)}$

Isso produz os parâmetros de uma distribuição normal aproximada. Para alta precisão, este Normal aproximado domina exceto nas caudas extremas. (Quando k < exp ( 5 ) , pode ser necessário aumentar um pouco o pdf Normal para garantir a dominação.)$f(u)$$k \lt \exp(5)$

Depois de realizar esse trabalho preliminar para qualquer valor dado de e estimar uma constante M > 1 (como descrito abaixo), obter uma variável aleatória é uma questão de:$k$$M \gt 1$

1. Desenhe um valor da distribuição normal dominante g ( u ) .$u$$g(u)$

2. Se ou se uma nova variável uniforme X exceder f ( u ) / ( M g $u \lt 1$$X$ , retorne à etapa 1.$f(u)/(M g(u))$

3. Defina .$x = u^{3/2}$

O número esperado de avaliações de devido às discrepâncias entre g e f é apenas ligeiramente maior do que 1. (Algumas avaliações adicionais irão ocorrer devido à rejeição de variates menos do que 1 , mas mesmo quando k é tão baixa quanto 2 a frequência de tais ocorrências é pequena.)$f$$g$$f$$1$$k$$2$

Este gráfico mostra os logaritmos de g e f em função de u para . Como os gráficos são muito próximos, precisamos inspecionar sua proporção para ver o que está acontecendo:$k=\exp(5)$

Isso exibe o log de razão de ; o fator M = exp ( 0,004 ) foi incluído para garantir que o logaritmo seja positivo em toda a parte principal da distribuição; isto é, para garantir M g ( u ) ≥ f ( u ), exceto possivelmente em regiões com probabilidade insignificante. Ao fazer M suficientemente grande, você pode garantir que M ⋅ $\log(\exp(0.004)g(u)/f(u))$$M = \exp(0.004)$$Mg(u) \ge f(u)$$M$$M \cdot g$ dominar $f$ em todas as caudas, exceto as mais extremas (que praticamente não têm chance de serem escolhidas em uma simulação). No entanto, quanto maior o , mais freqüentemente ocorrerão rejeições. À medida que k cresce, M pode ser escolhido muito próximo de 1 , o que implica praticamente nenhuma penalidade.$M$$k$$M$$1$

A similar approach works even for $k \gt \exp(2)$, but fairly large values of $M$ may be needed when $\exp(2) \lt k \lt \exp(5)$, because $f(u)$ is noticeably asymmetric. For instance, with $k = \exp(2)$, to get a reasonably accurate $g$ we need to set $M=1$:

A curva vermelha superior é o gráfico de enquanto a curva azul inferior é o gráfico de log ( f ( u ) ) . A amostragem de rejeição de f em relação a exp ( 1 ) g fará com que cerca de 2/3 de todos os sorteios sejam rejeitados, triplicando o esforço: ainda não é ruim. A cauda direita ( u > 10 ou x > 10 3 / ) estará sub-representada na amostragem de rejeição (porque $\log(\exp(1)g(u))$$\log(f(u))$$f$$\exp(1)g$$u \gt 10$$x \gt 10^{3/2} \sim 30$ não domina mais f ), mas essa cauda compreende menos que exp ( - 20 ) ∼ 10 - 9 da probabilidade total.$\exp(1)g$$f$$\exp(-20) \sim 10^{-9}$

To summarize, after an initial effort to compute the mode and evaluate the quadratic term of the power series of $f(u)$ around the mode--an effort that requires a few tens of function evaluations at most--you can use rejection sampling at an expected cost of between 1 and 3 (or so) evaluations per variate. The cost multiplier rapidly drops to 1 as $k = c d$ increases beyond 5.

Even when just one draw from $f$ is needed, this method is reasonable. It comes into its own when many independent draws are needed for the same value of $k$, for then the overhead of the initial calculations is amortized over many draws.

Addendum

@Cardinal has asked, quite reasonably, for support of some of the hand-waving analysis in the forgoing. In particular, why should the transformation $x = u^{3/2}$ make the distribution approximately Normal?

In light of the theory of Box-Cox transformations, it is natural to seek some power transformation of the form $x = u^\alpha$ (for a constant $\alpha$, hopefully not too different from unity) that will make a distribution "more" Normal. Recall that all Normal distributions are simply characterized: the logarithms of their pdfs are purely quadratic, with zero linear term and no higher order terms. Therefore we can take any pdf and compare it to a Normal distribution by expanding its logarithm as a power series around its (highest) peak. We seek a value of $\alpha$ that makes (at least) the third power vanish, at least approximately: that is the most we can reasonably hope that a single free coefficient will accomplish. Often this works well.

But how to get a handle on this particular distribution? Upon effecting the power transformation, its pdf is

Take its logarithm and use Stirling's asymptotic expansion of $\log(\Gamma)$:

(for small values of $c$, which is not constant). This works provided $\alpha$ is positive, which we will assume to be the case (for otherwise we cannot neglect the remainder of the expansion).

Compute its third derivative (which, when divided by $3!$, will be the coefficient of the third power of $u$ in the power series) and exploit the fact that at the peak, the first derivative must be zero. This simplifies the third derivative greatly, giving (approximately, because we are ignoring the derivative of $c$)

When $k$ is not too small, $u$ will indeed be large at the peak. Because $\alpha$ is positive, the dominant term in this expression is the $2\alpha$ power, which we can set to zero by making its coefficient vanish:

That's why $\alpha = 3/2$ works so well: with this choice, the coefficient of the cubic term around the peak behaves like $u^{-3}$, which is close to $\exp(-2 k)$. Once $k$ exceeds 10 or so, you can practically forget about it, and it's reasonably small even for $k$ down to 2. The higher powers, from the fourth on, play less and less of a role as $k$ gets large, because their coefficients grow proportionately smaller, too. Incidentally, the same calculations (based on the second derivative of $log(f(u))$ at its peak) show the standard deviation of this Normal approximation is slightly less than $\frac{2}{3}\exp(k/6)$, with the error proportional to $\exp(-k/2)$.

I like @whuber's answer very much; it's likely to be very efficient and has a beautiful analysis. But it requires some deep insight with respect to this particular distribution. For situations where you don't have that insight (so for different distributions), I also like the following approach which works for all distributions where the PDF is twice differentiable and that second derivative has finitely many roots. It requires quite a bit of work to set up, but then afterwards you have an engine that works for most distributions you can throw at it.

Basically, the idea is to use a piecewise linear upper bound to the PDF which you adapt as you are doing rejection sampling. At the same time you have a piecewise linear lower bound for the PDF which prevents you from having to evaluate the PDF too frequently. The upper and lower bounds are given by chords and tangents to the PDF graph. The initial division into intervals is such that on each interval, the PDF is either all concave or all convex; whenever you have to reject a point (x, y) you subdivide that interval at x. (You can also do an extra subdivision at x if you had to compute the PDF because the lower bound is really bad.) This makes the subdivisions occur especially frequently where the upper (and lower) bounds are bad, so you get a really good approximation of your PDF essentially for free. The details are a little tricky to get right, but I've tried to explain most of them in this series of blog posts - especially the last one.

Those posts don't discuss what to do if the PDF is unbounded either in domain or in values; I'd recommend the somewhat obvious solution of either doing a transformation that makes them finite (which would be hard to automate) or using a cutoff. I would choose the cutoff depending on the total number of points you expect to generate, say N, and choose the cutoff so that the removed part has less than $1 / (10 N)$ probability. (This is easy enough if you have a closed form for the CDF; otherwise it might also be tricky.)

This method is implemented in Maple as the default method for user-defined continuous distributions. (Full disclosure - I work for Maplesoft.)

I did an example run, generating 10^4 points for c = 2, d = 3, specifying [1, 100] as the initial range for the values:

There were 23 rejections (in red), 51 points "on probation" which were at the time in between the lower bound and the actual PDF, and 9949 points which were accepted after checking only linear inequalities. That's 74 evaluations of the PDF in total, or about one PDF evaluation per 135 points. The ratio should get better as you generate more points, since the approximation gets better and better (and conversely, if you generate only few points, the ratio is worse).

Você pode fazer isso executando numericamente o método de inversão, que diz que se você conectar variáveis ​​aleatórias uniformes (0,1) no CDF inverso, obterá um empate na distribuição. Incluí abaixo um código R que faz isso e, das poucas verificações que fiz, está funcionando bem, mas é um pouco desleixado e tenho certeza que você pode otimizá-lo.

Se você não está familiarizado com R, lgamma () é o log da função gama; integra () calcula uma integral 1-D definida; uniroot () calcula a raiz de uma função usando a bissecção 1-D.

# density. using the log-gamma gives a more numerically stable return for 
# the subsequent numerical integration (will not work without this trick)
f = function(x,c,d) exp( x*log(c) + (x-1)*log(d) - lgamma(x) )

# brute force calculation of the CDF, calculating the normalizing constant numerically
F = function(x,c,d) 
{
   g = function(x) f(x,c,d)
   return( integrate(g,1,x)$val/integrate(g,1,Inf)$val )
}

# Using bisection to find where the CDF equals p, to give the inverse CDF. This works 
# since the density given in the problem corresponds to a continuous CDF. 
F_1 = function(p,c,d) 
{
   Q = function(x) F(x,c,d)-p
   return( uniroot(Q, c(1+1e-10, 1e4))$root )
}

# plug uniform(0,1)'s into the inverse CDF. Testing for c=3, d=4. 
G = function(x) F_1(x,3,4)
z = sapply(runif(1000),G)

# simulated mean
mean(z)
[1] 13.10915

# exact mean
g = function(x) f(x,3,4)
nc = integrate(g,1,Inf)$val
h = function(x) f(x,3,4)*x/nc
integrate(h,1,Inf)$val
[1] 13.00002 

# simulated second moment
mean(z^2)
[1] 183.0266

# exact second moment
g = function(x) f(x,3,4)
nc = integrate(g,1,Inf)$val
h = function(x) f(x,3,4)*(x^2)/nc
integrate(h,1,Inf)$val
[1] 181.0003

# estimated density from the sample
plot(density(z))

# true density 
s = seq(1,25,length=1000)
plot(s, f(s,3,4), type="l", lwd=3)

A principal coisa arbitrária que faço aqui é assumir que $(1,10000)$é um suporte suficiente para a bissecção - fiquei preguiçoso com isso e pode haver uma maneira mais eficiente de escolher esse suporte. Para valores muito grandes, o cálculo numérico do CDF (por exemplo,$> 100000$) falhar, portanto o suporte deve estar abaixo disso. O CDF é efetivamente igual a 1 nesses pontos (a menos que$c, d$são muito grandes), então provavelmente poderia ser incluído algo que impediria o erro de cálculo do CDF para valores de entrada muito grandes.

Editar: quando$cd$é muito grande, ocorre um problema numérico com esse método. Como whuber aponta nos comentários, uma vez que isso ocorreu, a distribuição é essencialmente degenerada no seu modo, tornando-o um problema trivial de amostragem.