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

r=1, s=0.1

Below: r=3, s=1

r=3, s=1

Below: r=1, s=0

r=1, s=0

Below: r=100, s=5

r=100, s=5

• Luis Mendo: the distance to the origin follows a normal (Gaussian) distribution This is confusing, because any Gaussian distribution can produce negative numbers, and a distance cannot be negative. Your methods 1 and 2 (I haven't looked at what method 3 does) correspond to taking the absolute value of the Gaussian (and shifting the phase by 180 degrees, which is not significant)
• Jitse: @LuisMendo That's a valid point. The third method does basically the same. Consider the distance to be the displacement in the chosen angular direction.
• Nitrodon: @LuisMendo method 3 seems to be the same as methods 1 and 2, with the random angle generated as the direction of the vector formed by two iid normal variables.
• Andrey Tyukin: re "lack of a better term": in math (complex analysis etc.) that's usually called "annulus".

Answer 1

APL(Dyalog Unicode), ^{(削除) (削除ここまで)}38 bytes ^SBCS

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

{(×ばつ(.5*⍨ ×ばつ⍟?0)×ばつ1○しろまる○しろまる×ばつ?0)×ばつ ̄12○しろまる×ばつ○しろまる?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*⍨ ×ばつ⍟?0)×ばつ1○しろまる○しろまる×ばつ?0) draws a normally distributed value \$b \sim N(\alpha, \omega^2)\$:
?0 random number \$c \in (0,1)\$
1○しろまる○しろまる×ばつ?0: \$\sin(2\pi c)\$
?0 random number \$d \in (0,1)\$
.5*⍨ ×ばつ⍟?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○しろまる×ばつ○しろまる?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*⍨ ×ばつ⍟?0)×ばつ1○しろまる○しろまる×ばつ?0)×ばつ ̄12○しろまる×ばつ○しろまる?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')

enter image description here

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:

enter image description here

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

• user2554330: You could save a few bytes by using the R 4.1.0+ shorthand function definition notation: \(r,s)rnorm(1,r,s)*1i^runif(1,,4), though TIO doesn't seem to support it.
• Dominic van Essen: @user2554330 - Yes, but see comments here. TLDR: changing a single function into \ to save 7 bytes seems boring and I don't usually bother, since it's nice to have a TIO link.

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):
enter image description here

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 101円:

enter image description here

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:

enter image description here

r=3, s=1:

enter image description here

r=100, s=5:

enter image description here

Notes:

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

• att: The question asks for one point, not 1000; those 1000 points are generated for illustration. Random[] is a shorter way to access a uniform distribution. You can also save more bytes by using prefix/infix notation and ## instead of #,#2.

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) {}

• swpalmer: Try it online needs to update their JDK.
• Kevin Cruijssen: You should include the imports, so the var g=new Random(); should be var g=new java.util.Random();. But you can golf the var g=new java.util.Random();double a=g.nextGaussian(r,s),b=g.nextDouble(Math.PI*2); to double a=new java.util.Random().nextGaussian(r,s),b=Math.random()*Math.PI*2; for -8 bytes.
• swpalmer: @KevinCruijssen Thanks. updated. +2 bytes.. length still fits in Java's irritating signed byte 😂
• Kevin Cruijssen: Two more minor things to golf: (r,s)-> can be r->s-> for -1 byte and return new P(a*Math.cos(b),a*Math.sin(b)); can be return a*Math.cos(b)+","+a*Math.sin(b); for another -3 (with the Bifunction replaced with a currying Function: Function<Double, Function<Double, String>>)
• Kevin Cruijssen: And if you haven't seen them yet, tips for golfing Java and tips for golfing in 'all languages' might be interesting to read through. :)

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: Randomize. 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.

• emanresu A: τ is two bytes shorter than \tau. Sorry about the edit, my mistake\

Answer 12

05AB1E, 36 bytes

9°©D(ŸDÄ®>αÅ1* ̃Ω®/*+žq·®*ÝΩ®/DÅ3⁄4sÅ1⁄2‚*

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Å3⁄4sÅ1⁄2‚* # Calculate the resulting [x,y] pair using the two random
 # values:
D # Duplicate it
 Å3⁄4 # Pop and push its cosine
 s # Swap so the random value is at the top again
 Å1⁄2 # Pop and push its sine
 ‚ # Pair them together
 * # Multiply both to the earlier random value
 # (after which the result is output implicitly)