Simule uma distribuição uniforme em um disco

24

Eu estava tentando simular a injeção de pontos aleatórios dentro de um círculo, de forma que qualquer parte do círculo tivesse a mesma probabilidade de ter um defeito. Eu esperava que a contagem por área da distribuição resultante seguisse uma distribuição de Poisson se eu dividir o círculo em retângulos de área iguais.

Como requer apenas a colocação de pontos dentro de uma área circular, injetei duas distribuições aleatórias uniformes em coordenadas polares: (raio) e (ângulo polar).θRθ

Mas depois de fazer essa injeção, eu claramente recebo mais pontos no centro do círculo em comparação com a borda.

insira a descrição da imagem aqui

Qual seria a maneira correta de executar essa injeção no círculo, de modo que os pontos sejam distribuídos aleatoriamente no cirlce?

Jonjilla
fonte
Esta pergunta possui um analogia exata no fórum de Geometria: math.stackexchange.com/questions/87230/…
Aksakal

Respostas:

35

Você deseja que a proporção de pontos seja uniformemente proporcional à área e não à distância da origem. Como a área é proporcional à distância ao quadrado, gere raios aleatórios uniformes e tome suas raízes quadradas. Combine isso com um ângulo polar uniforme.

É rápido e simples de codificar, eficiente na execução (especialmente em uma plataforma paralela) e gera exatamente o número prescrito de pontos.

Exemplo

Este é um Rcódigo funcional para ilustrar o algoritmo.

n <- 1e4
rho <- sqrt(runif(n))
theta <- runif(n, 0, 2*pi)
x <- rho * cos(theta)
y <- rho * sin(theta)
plot(x, y, pch=19, cex=0.6, col="#00000020")

insira a descrição da imagem aqui

whuber
fonte
3

A amostragem por rejeição pode ser usada. Isso significa que podemos amostrar a partir da distribuição uniforme 2D e selecionar amostras que atendam às condições do disco.

Aqui está um exemplo.

x=runif(1e4,-1,1)
y=runif(1e4,-1,1)

d=data.frame(x=x,y=y)
disc_sample=d[d$x^2+d$y^2<1,]
plot(disc_sample)

insira a descrição da imagem aqui

Haitao Du
fonte
3
Essa é uma boa alternativa à abordagem adotada pelo OP. Simples e eficiente. No entanto, ele realmente não aborda a questão, que diz respeito a como modificar o método de coordenadas polares para produzir variáveis ​​uniformemente distribuídas. Por que podemos nos importar? Devido às implicações: depois de saber como gerar pontos uniformemente distribuídos nas coordenadas polares, você pode usar a amostragem por rejeição (e outros métodos familiares) nas coordenadas polares para amostrar de regiões que podem ser proibitivamente complicadas de serem amostradas em coordenadas cartesianas (pense em hipocicloides) , por exemplo).
whuber
1
π/4
@ whuber obrigado por me educar, comentando a minha resposta!
Haitao Du
3

Darei a você uma resposta n-dimensional geral que também funciona para casos bidimensionais, é claro. Em três dimensões, um análogo de um disco é o volume de uma bola sólida (esfera).

Existem duas abordagens que vou discutir. Um deles eu chamaria de "preciso" , e você obterá uma solução completa com ele em R. O segundo eu chamo de heurística , e é apenas a idéia, nenhuma solução completa é fornecida.

Solução "precisa"

Minha solução é baseada nos trabalhos de Marsaglia e Muller . Basicamente, isso acontece para que o vetor gaussiano normalizado de acordo com sua norma lhe desse pontos distribuídos uniformemente em uma hiperesfera d-dimensional:

insira a descrição da imagem aqui

d1/d

n <- 1e4
rho <- sqrt(runif(n))
# d - # of dimensions of hyperdisk
d = 2
r = matrix(rnorm(n*d),nrow=n,ncol=d)
x = r/rep(sqrt(rowSums(r^2))/rho,1)
plot(x[,1], x[,2], pch=19, cex=0.6, col="#00000020")

insira a descrição da imagem aqui

Aqui está um trecho de código para o caso 3d, ou seja, uma bola sólida:

library(scatterplot3d)
n <- 1e3
# d - # of dimensions of hyperdisk

d=3
rho <- (runif(n))^(1/d)
r = matrix(rnorm(n*d),nrow=n,ncol=d)
x = r/rep(sqrt(rowSums(r^2))/rho,1)

scatterplot3d(x[,1], x[,2], x[,3])

insira a descrição da imagem aqui

Abordagem heurística

Vn(R)=πn2Γ(n2+1)Rn
Rn

Eu=1dxEu2<R2

1d+2

Aksakal
fonte
@ Silverfish, você está certo, eu consertei o idioma
Aksakal
@Silverfish, ele é lento devido ao uso de Gaussian variates, mas poderia ser mais rápido do que simples amostragem rejeição em caso dimensional elevada, o que não é óbvio para muitos, mas é um assunto diferente
Aksakal
1/d,d
@ whuber, eu estava colando, corrigi um erro de digitação no poder do cubo. Se usarmos gaussiano, a amostragem por rejeição não será melhor; portanto, teremos que usar algo em forma de sino mais rápido que o gaussiano, você está certo
Aksakal
0

Aqui está uma solução alternativa em R:

n <- 1e4
## r <- seq(0, 1, by=1/1000)
r <- runif(n)
rho <- sample(r, size=n, replace=T, prob=r)
theta <- runif(n, 0, 2*pi)
x <- rho * cos(theta)
y <- rho * sin(theta)
plot(x, y, pch=19, cex=0.6, col="#00000020")

insira a descrição da imagem aqui

Q_Li
fonte
4
Você pode explicar esta resposta em inglês simples? Na verdade, não somos um site de ajuda para códigos, e as respostas somente com códigos devem ser desencorajadas.
gung - Restabelece Monica
5
0 01r <- seq(0, 1, by=1/10)
1
@whuber Obrigado por apontar isso. Na verdade, é minha principal idéia da solução. Minha abordagem foi gerar muitos círculos uniformes com raios variados e, para cada círculo, o número de pontos é proporcional ao comprimento de seu raio. Portanto, em um comprimento unitário de círculos com raios diferentes, o número de pontos é o mesmo. Para evitar a natureza discreta, poderíamos amostrar rde Uniform (0,1).
Q_Li 9/10/17