Cette méthode utilise Python et un algorithme de rejet.

Commencez par créer des variables uniformément réparties dans un carré. Puis rejetez tous les points qui se trouvent en dehors du demi-cercle qui vous intéresse. Ici, j'ai choisi le cas le plus simple, mais il peut être étendu à n'importe quel demi-cercle si vous ajoutez la rotation, la mise à l'échelle et la translation de ces points.

Le code est le suivant:

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()

Comme il utilise numpy pour toutes les étapes, c'est relativement rapide. Pour obtenir un nombre fixe de points, générez plus que ce dont vous avez besoin au départ et réduisez la taille du tableau. Depuis que j'ai commencé avec un carré parfait (aire = 4) et que j'ai seulement besoin d'un demi-cercle (aire = pi / 2), nous devons générer environ 2,6 fois plus de points que nécessaire à la fin.

Le La boucle for en python n'est pas très rapide, alors essayez de vous en tenir à numpy uniquement les fonctions et opérations.

Je ne sais pas ce que vous entendez par «uniformément».

Voici une approche pour générer des points distribués «uniformément» le long des axes x et y, mais pas très efficace

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()

Vous devez générer des points dans le rectangle englobant et supprimer les points qui ne sont pas dans le demi-cercle

# 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)

Vous devez générer des angles phi uniformément répartis et prendre le sqrt du rayon uniformément généré r (qui tient compte de ce que nous voulons échantillon uniformément sur la zone , voir l'explication ci-dessous), pour vous assurer que vous échantillonnez des points uniformément dans le demi-cercle.

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)

Explication

Pour générer des points uniformément répartis sur un (demi-) cercle, nous devons nous assurer que chaque aire ou segment infinitésimal est "frappé" avec la même probabilité. Nous pouvons simplement échantillonner phi à partir d'une distribution aléatoire uniforme [0, 1) , multiplié par np.pi (donc [0, pi) ), car tous les angles doivent avoir la même probabilité d'être échantillonnés. Mais si nous échantillonnons r à partir d'une distribution aléatoire uniforme dans [0, 1) , nous générons trop de points à petite radia, et pas assez à grande radia, puisque le area se développe comme r ** 2 . Pour tenir compte de ce fait, nous devons biaiser notre radia échantillonnée en conséquence, et dans ce cas, nous pouvons appliquer la polarisation en prenant simplement la racine carrée ( np.sqrt ), pour appliquer la pondération correcte à les valeurs de rayon échantillonnées, et prenez en compte la plus grande surface des anneaux extérieurs.

Une explication bien meilleure et plus approfondie se trouve ici: https://stackoverflow.com/a/50746409/1170207

Comparaison des performances des méthodes d'échantillonnage de rejet

Cette méthode étant essentiellement un échantillonnage par inversion méthode, nous comparons sa performance à un algorithme d'échantillonnage de rejet.

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()

Avec la méthode d'échantillonnage de rejet, nous ne pouvons pas non plus garantir que nous tirons un certain nombre de variables, nous devons donc répéter le processus jusqu'à ce que nous ayons . Cela ne peut pas être vectorisé aussi bien, sauf si nous acceptons d'échantillonner trop de valeurs et d'en supprimer d'autres.

Dans R:

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

Nous pouvons vérifier que cela fonctionne en évaluant une intégrale.

# 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

Nous évaluons maintenant numériquement l'intégrale de f .

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)