Random Point from a 2D Donut Distribution

Question

A donut distribution (for lack of a better term) is a random distribution of points in a 2-dimensional plane, forming a donut-like shape. The distribution is defined by two parameters: the radius r and spread s, in which the distance to the origin follows a normal (Gaussian) distribution around r, with a standard deviation s. The angular distribution is uniform in the range [0,2π).

The challenge

Given a radius r and spread s, your code should yield the Cartesian ((x,y)) coordinates of a single point chosen from this distribution.

Remarks

• Running your code multiple times with the same input should result in the specified distribution.
• Outputting polar coordinates is too trivial and not allowed.
• You can output Cartesian coordinates in any way allowed by the default I/O rules.
  • This includes complex values.

Valid approaches

Several algorithms can be used to yield the desired distribution, including but not limited to

1. Choose a from the uniform distribution [0,2π) and b from the normal distribution (r,s).
   Let x = b*cos(a) and y = b*sin(a).
2. Choose a from the uniform distribution [0,4) and b from the normal distribution (r,s).
   Let x+y*i = b*i^a.
3. Choose a,b,c all from the normal distribution (0,1).
   Let d = a+b*i and x+y*i = d/abs(d) * (c*s+r).

Example distributions (N=1000)

Below: r=1, s=0.1

Below: r=3, s=1

Below: r=1, s=0

Below: r=100, s=5

Answer 1

APL(Dyalog Unicode), 38 bytes ^SBCS

Takes r on the left and s on the right and returns a complex number.

{(⍺+⍵×(.5*⍨¯2×⍟?0)×1○○2×?0)ׯ12○2×○?0}

Try it on APLgolf!

Uses the first approach presented in the question. Dyalog doesn't have a builtin for sampling from a normal distribution, so this uses the Box–Muller transform to convert to random numbers from \$(0,1)\$ to a normally distributed value:

(⍺+⍵×(.5*⍨¯2×⍟?0)×1○○2×?0) draws a normally distributed value \$b \sim N(\alpha, \omega^2)\$:
?0 random number \$c \in (0,1)\$
1○○2×?0: \$\sin(2\pi c)\$
?0 random number \$d \in (0,1)\$
.5*⍨¯2×⍟?0: \$\sqrt{-2\ln{d}}\$
⍺+⍵× Scale from \$N(0, 1)\$ to \$N(\alpha, \omega^2)\$

?0 generates a random number \$a \in (0,1)\$
¯12○2×○?0: \$e^{i 2\pi a} = \sin(2\pi a)i + cos(2\pi a)\$

The product of these two values is the result.

Plotting code and images:

'InitCauseway' ⎕CY 'sharpplot'
InitCauseway ⍬
sp←⎕NEW Causeway.SharpPlot(700)
sp.SetTrellis(2 2)
sp.TrellisStyle←4

F ← {(⍺+⍵×(.5*⍨¯2×⍟?0)×1○○2×?0)ׯ12○2×○?0}

:For r s :In (1 0.1)(3 1)(1 0)(100 5)
    sp.NewCell
    sp.Heading←'r = ',(⍕r),'; s = ',⍕s
    sp.SetAxesScales(1)
    sp.DrawScatterPlot↓9 11∘.○{r F s}¨⍳1000
:EndFor

sp.SaveSvg(⊂'plot.svg')

Answer 2

MATL, 10 bytes

Xr*+Jr4*^*

Inputs are r, then s. Output is a complex number.

Try it online! Or see the plot for 1000 points at MATL Online! (it takes 10‒15 seconds).

How it works

Uses method 2 described in the challenge.

Xr   % Push random number with standard Gaussian distribution
*    % Implicit input: r. Multiply
+    % Implicit input: s. Add
J    % Push imaginary unit
r    % Push random number with stantard uniform distribution
4    % Push 4
*    % Multiply
^    % Power
*    % Multiply. Implicit output

Answer 3

Factor, 56 bytes

[ normal-random-float 2pi random 2dup cos * -rot sin * ]

Try it online!

Verifying correctness:

Explanation

                      ! 1 0.1
normal-random-float   ! 1.091729295255315
2pi                   ! 1.091729295255315 6.283185307179586
random                ! 1.091729295255315 4.669140230445313
2dup                  ! 1.091729295255315 4.669140230445313 1.091729295255315 4.669140230445313
cos                   ! 1.091729295255315 4.669140230445313 1.091729295255315 -0.04323526873134136
*                     ! 1.091729295255315 4.669140230445313 -0.04720120946224146
-rot                  ! -0.04720120946224146 1.091729295255315 4.669140230445313
sin                   ! -0.04720120946224146 1.091729295255315 -0.9990649185802335
*                     ! -0.04720120946224146 -1.090708439475907

Answer 4

R, 40 bytes

function(r,s)rnorm(1,r,s)*1i^runif(1,,4)

Try it online! or plot the results at rdrr.io

Answer 5

Wolfram Language (Mathematica), 42 bytes

Re[#+#2√Log[16/r^2]I^r]I^r&
r:=4Random[]

Try it online!

RandomVariate@NormalDistribution is costly (and, as noted by Ben Izd, doesn't work with stdev=0), so this uses Box-Muller to generate a normal distribution from two uniform ones.

Sample distributions (N=10000):

Answer 6

Julia 1.0, 26 bytes

r\s=(randn()s+r)im^4rand()

Try it online!

uses the second formula. output is a complex number. randn gives a random number from a normal distribution (0,1), and rand from a uniform distribtion in [0,1)

1000 points from 10\1:

Answer 7

Mathematica - 108 bytes

Following the first method, we could have:

fn=AngleVector/@RandomVariate[ProductDistribution[NormalDistribution[#,#2],UniformDistribution[{0,2Pi}]],1000]&;

then visualize it by:

visualize = 
  Graphics[{PointSize[.01], Point[fn[#, #2]]}, Frame -> True] &;

r=1, s=0.1:

r=3, s=1:

r=100, s=5:

Notes:

• Variance: 0 in NormalDistribution is not supported (could be hacked by having a small number)

Answer 8

Python 2, 58 bytes

Might not be the shortest, but here is the base-case for python i guess.

Outputs a complex number:

lambda r,s:1j**uniform(0,4)*gauss(r,s)
from random import*

Try it online!

Answer 9

R, 45 bytes

function(r,s)exp(runif(1)*2i*pi)*rnorm(1,r,s)

Try it online!

Uses the first method in the description, of course using the neat fact that \$e^{i\theta}=\cos(\theta)+i\sin(\theta)\$. Longer than Dominic van Essen's answer by 5 bytes, though.

Answer 10

Java 17, 127 bytes

(r,s)->{double a=new java.util.Random().nextGaussian(r,s),b=Math.random()*Math.PI*2;return new P(a*Math.cos(b),a*Math.sin(b));}

This is a BiFunction<Double, Double, P> where P is a record P(double x, double y) {}

Answer 11

Desmos, 80 70 68 bytes

b=random(1,t)τ
f(r,s,t)=(normaldist(r,s).random(1,t)(cosb,sinb))[1]

Takes an extra argument t as the seed, which is the only way to re-use the function to get different samples without pressing the randomize button.

Try it on Desmos!

-10 bytes thanks to Aiden Chow

-2 bytes thanks to emanrescu A (\tau to τ)

All functions in Desmos are pure, so they can't return a different value when evaluated at different times, even in a list comprehension. This causes an issue with the CGSE policy of functions being re-usable.

There's a randomize button to re-seed all of the random seed-dependent function calls: . This doesn't vibe with me because it requires user interaction to re-seed, but it would allow the following 53-byte submission:

b=random()τ
f(r,s)=(cosb,sinb)normaldist(r,s).random

In this submission, I opted to take the random seed as an extra argument, which is a common design decision in Desmos if a program needs to avoid user action when re-seeding. This is the only way to get different outputs from random functions without the user pressing the randomize button.

Answer 12

05AB1E, 36 bytes

9°©D(ŸDÄ®>αÅ1*˜Ω®/*+žq·®*ÝΩ®/DžsŽ‚*

Inputs in the order \$r,s\$.

Try it online. (With the 9 replaced with 3 so it won't time out.)

Explanation:

Uses the first method described in the challenge description.

However, 05AB1E lacks a Gaussian distribution random builtin, as well as a builtin to get a random decimal number given a range. Both of those are therefore done manually.

9°©D(ŸDÄ®>αÅ1*˜Ω®/ # Push a random value with Gaussian distribution within the
                   # range [-1,1]:
9°                 #  Push 1,000,000,000 (10**9)
  ©                #  Store it in variable `®` (without popping)
   D               #  Duplicate it
    (              #  Negate the copy
     Ÿ             #  Pop both and push an integer-list list in the range
                   #  [-1000000000,1000000000]
      D            #  Duplicate this list
       Ä           #  Get the absolute value of each in the copy
        ®>α        #  Get the absolute difference between each and `®`+1
           Å1      #  Map each inner value to a list of that many 1s
             *     #  Multiply each to the values at the same positions in the
                   #  remaining list
              ˜    #  Flatten this list of lists
               Ω   #  Pop and push a random integer
                ®/ #  Divide it by `®`
*+                 # Use the inputs to transform it into s*random+r:
*                  #  Multiply it to the (implicit) input `s`
 +                 #  Add the (implicit) input `r`
žq·®*ÝΩ®/          # Push a random value with uniform distribution within the
                   # range [0,2π):
žq                 #  Push builtin PI: 3.141592653589793
  ·                #  Double it to get tau: 6.283185307179586
   ®*              #  Multiply it by `®`
     Ý             #  Pop and push an integer list in the range [0,2π®]
      Ω            #  Pop and push a random integer from this list
       ®/          #  Divide it by `®`
DžsŽ‚*           # Calculate the resulting [x,y] pair using the two random
                   # values:
D                  #  Duplicate it
 ž                #  Pop and push its cosine
   s               #  Swap so the random value is at the top again
    Ž             #  Pop and push its sine
      ‚            #  Pair them together
       *           #  Multiply both to the earlier random value
                   # (after which the result is output implicitly)