Gere pontos com eficiência entre o círculo e o quadrado da unidade

Eu gostaria de gerar amostras da região azul definida aqui:

A solução ingênua é usar a amostragem por rejeição no quadrado da unidade, mas isso fornece apenas uma eficiência de (~ 21,4%).$1-\pi/4$

Existe alguma maneira de provar com mais eficiência?

Dois milhões de pontos por segundo serão suficientes?

A distribuição é simétrica: só precisamos calcular a distribuição para um oitavo do círculo completo e depois copiá-la em torno dos outros octantes. Nas coordenadas polares , a distribuição cumulativa do ângulo Θ para a localização aleatória ( X , Y ) no valor θ é dada pela área entre o triângulo ( 0 , 0 ) , ( 1 , 0 ) , ( 1 , tan θ ) e o arco do círculo que se estende de $(r,\theta)$$\Theta$$(X,Y)$$\theta$$(0,0), (1,0), (1,\tan\theta)$ a ( cos θ , sin θ ) . É assim proporcional a$(1,0)$$(\cos\theta,\sin\theta)$

$$F_\Theta(\theta) = \Pr(\Theta \le \theta) \propto \frac{1}{2}\tan(\theta) - \frac{\theta}{2},$$

de onde sua densidade é

$$f_\Theta(\theta) = \frac{d}{d\theta} F_\Theta(\theta) \propto \tan^2(\theta).$$

Podemos amostrar a partir dessa densidade usando, digamos, um método de rejeição (que tem eficiência ).$8/\pi-2 \approx 54.6479\%$

A densidade condicional da coordenada radial é proporcional à entre e . Isso pode ser amostrado com uma fácil inversão do CDF.r d r r = 1 r = seg $R$$rdr$$r=1$$r=\sec\theta$

Se amostras independentes , a conversão de volta para coordenadas cartesianas amostra desse octante. Como as amostras são independentes, a troca aleatória de coordenadas produz uma amostra aleatória independente do primeiro quadrante, conforme desejado. (Os swaps aleatórios requerem a geração de apenas uma única variável binomial para determinar quantas das realizações a serem trocadas.)( x i , y i$(r_i,\theta_i)$$(x_i,y_i)$

Cada realização de requer, em média, uma variável uniforme (para ) mais vezes duas variáveis ​​uniformes (para ) e uma pequena quantidade de cálculo (rápido). São variáveis ​​por ponto (que, é claro, tem duas coordenadas). Detalhes completos estão no exemplo de código abaixo. Esse número representa 10.000 dos mais de meio milhão de pontos gerados.R 1 / ( 8 π - 2 ) Θ 4 / ( π - 4 ) ≈ $(X,Y)$$R$$1/(8\pi-2)$$\Theta$$4/(\pi-4) \approx 4.66$

Aqui está o Rcódigo que produziu esta simulação e cronometrou-a.

n.sim <- 1e6
x.time <- system.time({
  # Generate trial angles `theta`
  theta <- sqrt(runif(n.sim)) * pi/4
  # Rejection step.
  theta <- theta[runif(n.sim) * 4 * theta <= pi * tan(theta)^2]
  # Generate radial coordinates `r`.
  n <- length(theta)
  r <- sqrt(1 + runif(n) * tan(theta)^2)
  # Convert to Cartesian coordinates.
  # (The products will generate a full circle)
  x <- r * cos(theta) #* c(1,1,-1,-1)
  y <- r * sin(theta) #* c(1,-1,1,-1)
  # Swap approximately half the coordinates.
  k <- rbinom(1, n, 1/2)
  if (k > 0) {
    z <- y[1:k]
    y[1:k] <- x[1:k]
    x[1:k] <- z
  }
})
message(signif(x.time[3] * 1e6/n, 2), " seconds per million points.")
#
# Plot the result to confirm.
#
plot(c(0,1), c(0,1), type="n", bty="n", asp=1, xlab="x", ylab="y")
rect(-1, -1, 1, 1, col="White", border="#00000040")
m <- sample.int(n, min(n, 1e4))
points(x[m],y[m], pch=19, cex=1/2, col="#0000e010")

Proponho a solução a seguir, que deve ser mais simples, mais eficiente e / ou computacionalmente mais barata do que outras almas de @cardinal, @whuber e @ stephan-kolassa até agora.

Envolve as seguintes etapas simples:

1) Desenhe duas amostras uniformes padrão:

$$u_1 \sim Unif(0,1)\\ u_2 \sim Unif(0,1).$$

2a) Aplique a seguinte transformação de cisalhamento ao ponto (os pontos no triângulo inferior direito são refletidos no triângulo superior esquerdo e estarão "des- refletido "em 2b): $\min\{u_1,u_2\}, \max\{u_1,u_2\}$

$$\begin{bmatrix} x\\y \end{bmatrix} = \begin{bmatrix} 1\\1 \end{bmatrix} + \begin{bmatrix} \frac{\sqrt{2}}{2} & -1\\ \frac{\sqrt{2}}{2} - 1 & 0\\ \end{bmatrix} \, \begin{bmatrix} \min\{u_1,u_2\}\\ \max\{u_1,u_2\}\\ \end{bmatrix}.$$

2b) Troque e se .$x$$y$$u_1 > u_2$

3) Rejeite a amostra se dentro do círculo de unidades (a aceitação deve ser em torno de 72%), ou seja:

$$x^2 + y^2 < 1.$$

A intuição por trás desse algoritmo é mostrada na figura.

As etapas 2a e 2b podem ser mescladas em uma única etapa:

2) Aplique a transformação de cisalhamento e troque

$$x = 1 + \frac{\sqrt{2}}{2} \min(u_1, u_2) - u_2\\ y = 1 + \frac{\sqrt{2}}{2} \min(u_1, u_2) - u_1$$

O código a seguir implementa o algoritmo acima (e o testa usando o código do @ whuber).

n.sim <- 1e6
x.time <- system.time({
    # Draw two standard uniform samples
    u_1 <- runif(n.sim)
    u_2 <- runif(n.sim)
    # Apply shear transformation and swap
    tmp <- 1 + sqrt(2)/2 * pmin(u_1, u_2)
    x <- tmp - u_2
    y <- tmp - u_1
    # Reject if inside circle
    accept <- x^2 + y^2 > 1
    x <- x[accept]
    y <- y[accept]
    n <- length(x)
})
message(signif(x.time[3] * 1e6/n, 2), " seconds per million points.")
#
# Plot the result to confirm.
#
plot(c(0,1), c(0,1), type="n", bty="n", asp=1, xlab="x", ylab="y")
rect(-1, -1, 1, 1, col="White", border="#00000040")
m <- sample.int(n, min(n, 1e4))
points(x[m],y[m], pch=19, cex=1/2, col="#0000e010")

Alguns testes rápidos produzem os seguintes resultados.

Algoritmo /stats//a/258349 . Melhor de 3: 0,33 segundos por milhão de pontos.

Este algoritmo. Melhor de 3: 0,18 segundos por milhão de pontos.

Bem, com mais eficiência, isso pode ser feito, mas espero que você não esteja procurando mais rapidamente .

A idéia seria provar primeiro um valor , com uma densidade proporcional ao comprimento da fatia azul vertical acima de cada valor :$x$$x$

$$f(x) = 1-\sqrt{1-x^2}.$$

A Wolfram ajuda você a integrar isso :

$$\int_0^x f(y)dy = -\frac{1}{2}x\sqrt{1-x^2}+x-\frac{1}{2}\arcsin x.$$

Portanto, a função de distribuição cumulativa seria essa expressão, dimensionada para integrar a 1 (isto é, dividida por ).∫ 1 0 f ( y ) d $F$$\int_0^1 f(y)dy$

Agora, para gerar seu valor , escolha um número aleatório , distribuído uniformemente entre e . Então encontre tal que . Ou seja, precisamos inverter o CDF ( amostragem por transformação inversa ). Isso pode ser feito, mas não é fácil. Nem rápido.t 0 1 x F ( x ) = $x$$t$$0$$1$$x$$F(x)=t$

Finalmente, dado , escolha um aleatório distribuído uniformemente entre e .y $x$$y$$\sqrt{1-x^2}$$1$

Abaixo está o código R. Observe que estou pré-avaliando o CDF em uma grade de valores , e mesmo assim isso leva alguns minutos.$x$

Provavelmente, você pode acelerar bastante a inversão do CDF se investir algum pensamento. Então, novamente, o pensamento dói. Eu, pessoalmente, iria para a amostragem de rejeição, o que é mais rápido e muito menos propenso erro de, a menos que eu tinha muito boas razões para não.

epsilon <- 1e-6
xx <- seq(0,1,by=epsilon)
x.cdf <- function(x) x-(x*sqrt(1-x^2)+asin(x))/2
xx.cdf <- x.cdf(xx)/x.cdf(1)

nn <- 1e4
rr <- matrix(nrow=nn,ncol=2)
set.seed(1)
pb <- winProgressBar(max=nn)
for ( ii in 1:nn ) {
    setWinProgressBar(pb,ii,paste(ii,"of",nn))
    x <- max(xx[xx.cdf<runif(1)])
    y <- runif(1,sqrt(1-x^2),1)
    rr[ii,] <- c(x,y)
}
close(pb)

plot(rr,pch=19,cex=.3,xlab="",ylab="")