restart:
The Triangle
First interpretation : a statistical point of view based on the Dirichlet distribution
with(Statistics):# sample a triple (x, y, z) according to a Dir(1, 1, 1) distribution
N := 10000:
u := Sample(GammaDistribution(1, 1), N):
v := Sample(GammaDistribution(1, 1), N):
w := Sample(GammaDistribution(1, 1), N):
s := u +~ v +~ w:
x := u /~ s:
y := v /~ s:
z := w /~ s:
M := < x^+ | y^+ | z^+ >:
pxyz := PLOT3D(POINTS(convert(M, listlist))):
pxy := PLOT(POINTS(convert(M[..,[1,2]], listlist)), SYMBOL(POINT, 1)):
pyz := PLOT(POINTS(convert(M[..,[2,3]], listlist)), SYMBOL(POINT, 1)):
pzx := PLOT(POINTS(convert(M[..,[3,1]], listlist)), SYMBOL(POINT, 1)):# DocumentTools:-Tabulate([[pxyz, pxy],[pyz, pzx]]):Histogram(x):
Second interpretation : a geometrical point of view
The nodes of the small triangles are regularly spreaded (not uniformly in a statistical sense)
So are their barycenters.
At the very end one can also uniformly sample a point within each sub triangle to obtain a quasi uniform distribution restart:
T := < <0, 1, 0, 0> | <0, 0, 1, 0> >;Hom := (a, b, M) -> < a *~ M[..,1] | b *~ M[..,2] >:Trans := (a, b, M) -> < a +~ M[..,1] | b +~ M[..,2] >:N := 20:
PLOT(CURVES(T), seq(seq(CURVES(Trans(k/N, l/N, Hom(1/N, 1/N, T))), k=0..N-l), l=0..N)):AllNodes := {seq(seq(op(convert(Trans(k/N, l/N, Hom(1/N, 1/N, T)), listlist)), k=0..N-l), l=0..N)}:
AllNodes := map(evalf, convert(AllNodes, list)):
AllNodes := convert(AllNodes, Matrix):Statistics:-Histogram(AllNodes[..,1], binbounds=[seq((k-1/2)/N, k=0..N+2)], frequencyscale=absolute):
The Disc
First interpretation : geometrical points of view
restart:
# compare this one
radii := [seq(0.1..1,0.1)]:
angles := [seq(0..evalf(2*Pi),evalf(Pi/10))]:
pRhoTheta := plots:-display(
seq(plottools:-circle([0,0],k), k in radii),
seq(plottools:-line([0,0],[cos(k), sin(k)]), k in angles),
seq(seq(plottools:-point([u*cos(v), u*sin(v)], color=red, symbol=solidcircle, symbolsize=15), u in radii), v in angles)
):
# to this one
x := [seq(-1..1,0.1)]:
y := x:
xy := [seq(seq([u,v], u in x), v in y)]:
xy := map(u -> if u[1]^2+u[2]^2 <=1 then u end if, xy):
pXY := PLOT(POINTS(xy, COLOR(RGB, 1, 0, 0), SYMBOL(_SOLIDCIRCLE, 15))):
# DocumentTools:-Tabulate([pRhoTheta, pXY]):# Wich of these two plots give the better representation of the naive idea of uniformity ?
# If you say it is the right one, that means you implicitly consider that "uniformity" (better to say here
# "regularity") is associated to an euclidian metric.
# But, if you do not care about the metric, the left picture is as regular than the right one is.
# So : how can we judge of the regularity (next "uniformity") od a collection of points ?
# Based on the idea used on the right figure one can also build a set of N bi-uniformly distributed points
# in the square [-1, 1]x[-1, 1] square and keep only those inside the disc.
# This is a rejection-acceptation technique (basically the one vv used for sampling the triangle [but here
#without any mirroring to have a null rejection ratio as vv did])
Second interpretation : a statistical point of view
restart:
with(Statistics):# We sample the disc in polar (u, v) coordinates (u=radius, v=angle)
# The density must be constant in each annulus of thickness t and mean radius R
# Because the surface of such annulus is 2*Pi*R*t, it is clear that the sampling
# of u must be done from a triangular distribution.
N := 10000:
u := Sample(TriangularDistribution(0, 1, 1), N):
v := Sample(UniformDistribution(0, 2*Pi), N):
M0 := < u^+ | v^+ >:
# prepare a cartesian plot
M := < M0[..,1] *~ cos~(M0[..,2]) | M0[..,1] *~ sin~(M0[..,2]) >: plots:-display(plottools:-circle([0,0],1), PLOT(POINTS(convert(M, listlist)), SYMBOL(POINT, 1))):
The Ellipse
A statistical point of view
Idea : start from the disc case and map the result onto an ellipserestart:
with(Statistics):
N := 10000:
u := Sample(TriangularDistribution(0, 1, 1), N):
v := Sample(UniformDistribution(0, 2*Pi), N):
MC := < u^+ | v^+ >:
A := 3 ; B := 1 ;
ME := < A *~ MC[..,1] *~ cos~(MC[..,2]) | B *~ MC[..,1] *~ sin~(MC[..,2]) >: plots:-display(plottools:-ellipse([0,0], A, B), PLOT(POINTS(convert(ME, listlist)), SYMBOL(POINT, 0,3), SCALING(CONSTRAINED))):It seems uniform but it's not, at least in some reasonable sense.As a matter of fact, consider a point P in the open ellipse and a neighbourhood D(P) of P sufficiently large for it contains at least a few points Q1, Q2, ...
It is obvious that the distance between P and each of these Qn depends on the orientation of the segment (P, Qn) ; here
ratio of this distance, when P and Qn are aligned along the X axis, to the distance when they are aligned along the Y
axis is roughly A/B
So the distance between 2 arbitrary choosen points is not isotropic
A uniform and isotropic distribution : the acceptation/rejection algorithm
Sample the circumscribed square [-A, A]x[-A,A] ans reject points outside the ellipse
(efficiency (4*A*B) / (Pi*A*A) = (4*B) / (Pi*A) ; here about 0.4
Nota : sampling the rectangle [-A, A]x[-B, B] would give a better efficiency but the scaling effect would still exist, leading
to an anisotropic distribution
restart:
with(Statistics):
A := 3 ; B := 1 ;
EfficiencyRatio := evalf ((4*A*B) / (Pi*A*A));
N := round(10000 / EfficiencyRatio):
u := Sample(UniformDistribution(-A, A), N):
v := Sample(UniformDistribution(-A, A), N):
MC := < u^+ | v^+ >:
# A point M inside the ellipse lust verufy d(M,F)+d(M,F') <= 2*A
C := evalf(sqrt(A^2-B^2));
F1M := sqrt~((MC[..,1] -~ C)^~2 +~ MC[..,2]^~2):
F2M := sqrt~((MC[..,1] +~ C)^~2 +~ MC[..,2]^~2):
F1MF2 := convert(F1M +~ F2M, list):
PointsIn := zip((u,v) -> if u < 2*A then v end if, F1MF2, [seq(1..N)]):
MCin := MC[PointsIn, ..]:
plots:-display(
plottools:-ellipse([0,0], 3, 1, color=red),
PLOT(POINTS(convert(MC, listlist), SYMBOL(POINT, 1)), SCALING(CONSTRAINED)),
PLOT(POINTS(convert(MCin, listlist), SYMBOL(POINT, 1)), SCALING(CONSTRAINED), COLOR(RGB, 1, 0, 0))
):Last but not least
The more efficient way to make a set of points "uniformly"and isotropically distributed within an ellipse
is based on tesselations.
Suppose you have created a "uniform" tesselation of the ellipse made only of equilateral triangles (such a function exists in
Matlab : distmesh2D) : then the barycenters of each of these triangles form a set of point with the desired propreties ...... provided that we interpret uniformity as I explain to Markiyan (discrepancy) in a previous answer