Как создать набор случайных точек в полукруге равномерно в Python или R?

Question

Я пытаюсь сгенерировать набор точек в полукруге равномерно.

size = 1000
t  = np.random.random(size)*np.pi*2
u  = np.random.random(size) + np.random.random(size)
r = np.where(u > 1, 2 - u, u)
x = r*cos(t)
y = r*sin(t)
coor = (x,y)
for idx, value in enumerate(y):
    if value<0:
        x[idx]=-3
        y[idx]=-3
f, ax = plt.subplots(figsize = (3,3))
plt.scatter(x, y)

этот фрагмент кода имеет 2 ошибки.

1. На графике есть много (-3, -3) точек, которые нужно удалить
2. это кажется неэффективным.

рисунок, показанный ниже, НЕ является равномерным, поскольку точки больше центра, чем другие.

другой график, показанный ниже, можно рассматривать как равномерно.

любая идея по исправлению ошибок будет принята.

Answer 1

Вы должны сгенерировать равномерно распределенные углы phi и взять sqrt равномерно сгенерированного радиуса r (который учитывает, что мы хотим произвести равномерную выборку в области , см. Пояснение ниже ), чтобы обеспечить равномерную выборку точек на полукруге.

import numpy as np
import matplotlib.pyplot as plt

# sample
size = 10000
R = 1
phi = np.random.random(size=size) * np.pi
r = np.sqrt(np.random.random(size=size)) * R

# transform
x = r * np.cos(phi)
y = r * np.sin(phi)

# plot
f = plt.figure(figsize=(12,12))
a = f.add_subplot(111)
a.scatter(x, y, marker='.')
a.set_aspect('equal')
plt.show()

Объяснение

Чтобы сформировать равномерно распределенные точки на (полукруге) окружности, мы должны убедиться, что каждая бесконечно малая область или сегмент "поражены" с одинаковой вероятностью. Мы можем просто выбрать phi из равномерного случайного распределения [0, 1), умноженного на np.pi (т. Е. [0, pi)), поскольку все углы должны иметь одинаковую вероятность для выборки. Но если мы выберем r из равномерного случайного распределения в [0, 1), мы создадим слишком много точек при малых радиусах и недостаточно при больших радиусах, поскольку область увеличивается как r**2. Чтобы принять этот факт во внимание, мы должны соответствующим образом сместить наши выборочные радиусы, и в этом случае мы можем применить смещение, просто взяв квадратный корень (np.sqrt), чтобы применить правильное взвешивание к выборочным значениям радиуса, и принять с учетом большей площади наружных колец.

Здесь гораздо лучше и более подробное объяснение: https://stackoverflow.com/a/50746409/1170207

Сравнение производительности с методами выборки отбраковки

Поскольку этот метод в основном является методом инверсионной выборки, мы сравниваем его производительность по алгоритму выборки отклонения.

import numpy as np
x, y = np.random.random(size=(2,10000))
%timeit r, phi = np.sqrt(x), y
# 19 µs ± 33.4 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
%timeit m = x**2 + y**2 <= 1; xx, yy = x[m], y[m]
# 81.5 µs ± 271 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

С помощью метода отбраковки выборки мы также не можем гарантировать, что мы нарисуем определенное количество вариантов, поэтому мы должны повторять процесс до тех пор, пока у нас не получится. Это не может быть векторизовано так хорошо, если мы не примем выборку слишком большого количества значений и не отбросим дополнительные.

Answer 2

В R:

runif_in_semicircle <- function(n, radius=1){
  theta <- runif(n, 0, pi)
  r <- radius * sqrt(runif(n))
  cbind(r*cos(theta), r*sin(theta))
}

sims <- runif_in_semicircle(1000)
plot(sims[,1], sims[,2], asp=1, pch=19)

Мы можем проверить это, оценивая интеграл.

# integrand example
f <- function(x) x[1]^2 + exp(x[2])

set.seed(666)
sims <- runif_in_semicircle(10000)
fsims <- apply(sims, 1, f)
mean(fsims)*pi/2 # approximates the integral of f over the half-disk
# 2.890905

Сейчасмы численно оцениваем интеграл от f.

library(SphericalCubature)
adaptIntegrateBallPolar(f, n=2, lowerLimit = 0, upperLimit = pi)
# $integral
# [1] 2.880598

Answer 3

Вы должны сгенерировать точки в окружающем прямоугольнике и удалить точки, которые не находятся в полукруге

# generate about n points in a half-circle with
# given center (x, y) and given radius, and y>0

points <- function(x, y, radius, n) {
    n2 = n * 4 / pi # each point has pi/4 probability to survive

    # make [-1, 1] * [-1, 1] square
    xs = runif(n2, -1, 1)
    ys = runif(n2, 0, 1)  # or just runif(n2)
    points = cbind(xs, ys)

    # keep only points in circle with center (0,0) and radius 1 with y>0
    ind = (xs**2 + ys**2 <= 1) # the condition ys>=0 is obeyed already
    points = points[ind,]

    # move/stretch to given center and radius
    points = points * radius
    points[,1] = points[,1] + x
    points[,2] = points[,2] + y
}

# generate about 1000 points with center(1,1) and radius 3
points = f(1, 1, 3, 1000)

# plot them, making them smaller for better visibility
plot(points, cex=0.3)

Answer 4

Не уверен, что вы имеете в виду «единообразно».

Вот один из подходов для генерации «равномерно» распределенных точек вдоль осей x и y, но не очень эффективный

import numpy as np
import matplotlib.pyplot as plt

x = np.random.random(size) * 2 - 1
y = np.random.random(size)
r = np.sqrt(x**2 + y**2)
x[r > 1] = -3
y[r > 1] = -3
plt.plot(x, y, 'o')
plt.show()

Answer 5

Этот метод использует Python и алгоритм отклонения.

Сначала создайте равномерно распределенные переменные в квадрате. Затем отбросьте все точки, которые лежат за пределами интересующего вас полукруга. Здесь я выбрал самый простой случай, но он может быть расширен до любого полукруга, если вы добавите вращение, масштабирование и перемещение этих точек.

Код выглядит следующим образом:

import numpy as np
import matplotlib.pyplot as plt

n = 10000
r = 3  # radius

uniform_square = np.random.uniform(-r,r,(int(2.65*n),2))
radius = np.sqrt(uniform_square[:,0]**2 + uniform_square[:,1]**2)

uniform_circle = uniform_square[radius<=r,:]

uniform_half_circle = uniform_circle[uniform_circle[:,0]>=0,:]

final_points = uniform_half_circle[0:n,:]

fig, axs = plt.subplots(1, 1)
axs.scatter(final_points[:,0], final_points[:,1], marker='.')
axs.set_aspect('equal', 'box')
plt.show()

Поскольку он использует numpy для всех шагов, это относительно быстро. Чтобы получить фиксированное количество точек, сначала сгенерируйте больше, чем вам нужно, и уменьшите массив. Так как я начал с идеального квадрата (площадь = 4) и мне нужен только полукруг (область = пи / 2), нам нужно сгенерировать примерно в 2,6 раза больше точек, чем нам нужно в конце.

Цикл for в python не очень быстрый, поэтому старайтесь придерживаться numpy только функций и операций.

Answer 6

Вы можете попробовать следующее в R:

size <- 1000
maxRad <- 1 #maximum radius of the half-circle
r <- runif(1000,0,maxRad) #generate random radius
phi <- runif(size,0,pi) #generate angle for polarcoordinates (between 0 and pi since you want a halfcircle)
x <- r*cos(phi) #polarcoordinates
y <- r*sin(phi)
plot(x,y)

Вы можете поместить его в функцию

halfCircle <- function(size, maxRad) {
 r <- runif(1000,0,maxRad)
 phi <- runif(size,0,1)*pi
 x <- r*cos(phi)
 y <- r*sin(phi)
 plot(x,y)
}

и проверьте, дает ли он вам "более приемлемые случайные" результаты.