How to choose random float in a square weighted against arbitrary zones?

220 views Asked by At

I have a square, which goes from -1 to 1 in x and y.

Choosing a random point in this square is pretty easy:

Random r = new Random();
float x = (float)Math.Round(r.NextDouble() * 2 - 1, 4);
float y = (float)Math.Round(r.NextDouble() * 2 - 1, 4);

This gives me any point, with equal probability, in my square.

enter image description here

It woud also be pretty easy to just remove a section of the square from the possibilities

Random r = new Random();
float x = (float)Math.Round(r.NextDouble() * 1.5 - 1, 4);
float y = (float)Math.Round(r.NextDouble() * 2 - 1, 4);

enter image description here

But what I'm really struggling to do, is to weight the random towards a certain zone. Specifically, I would like the section highlighted here to be more likely, and everything else (except the red section, which is still off-limits) should have a probability lower depending on the distance from the highligthed line. The furthest point should have 0 chance, and the rest an existing chance which is higher when closer to the line, with points exactly on my line (since I round them to a specific decimal, there are points with are on the line) having the best odds.

enter image description here

Sorry for the ugly pictures. This is the best i could do in paint to show my thoughts.

The "most likely" area is an empty diamond (just the that with the vertices (-1, 0), (0, -0.5), (1, 0), (0, 0.5), with of course the red area override the weighting because it's off limits. The red area is anything with x > 0.5

Does anyone knows how to do this? I'm working in C# but honestly an algorithm in any non-esoteric language would do the trick. I'm completely lost as to how to proceed.

A commenter noted that adding the off-limits zone to the algorithm is an added difficulty with no real use.

You can assume that I'll take care of the off-limit section by myself after running the weighting algorithm. Since it's just 25% of the area, most of the times it wouldn't even make a difference performance-wise if I just made this:

while (x > 0.5)
{
    runAlgorithmAgain();
}

So you can safely ignore that part for answers.

1

There are 1 answers

0
Severin Pappadeux On BEST ANSWER

Ok, here my thoughts on this matter. I would like to propose algorithm which, with some rejections, might solve your problem. Note, due to need of acceptance-rejection, it might be slower than you expected it to be.

We sample in single quadrant (say, lower left one), then use reflection to put point into any other quadrant, and then reject red zone points.

Basically, sampling in quadrant is two-step process. First, we sample first position on the border line. As soon as we got position on the line, we sample from distribution which is bell-like shape (Gaussian or Laplace for example), and move point in the orthogonal to the border line direction.

Code compiles, but completely untested, so please check everything startign with the numbers

using System;

namespace diamond
{
    class Program
    {
        public const double SQRT_5 = 2.2360679774997896964091736687313;

        public static double gaussian((double mu, double sigma) N, Random rng) {
            var phi = 2.0 * Math.PI * rng.NextDouble();
            var r   = Math.Sqrt( -2.0 * Math.Log(1.0 - rng.NextDouble()) );
            return N.mu + N.sigma * r * Math.Sin(phi);
        }

        public static double laplace((double mu, double sigma) L, Random rng) {
            var v = - L.sigma * Math.Log(1.0 - rng.NextDouble());
            return L.mu + ((rng.NextDouble() < 0.5) ? v : -v );
        }

        public static double sample_length(double lmax, Random rng) {
            return lmax * rng.NextDouble();
        }

        public static (double, double) move_point((double x, double y) pos, (double wx, double wy) dir, double l) {
            return (pos.x + dir.wx * l, pos.y + dir.wy * l);
        }

        public static (double, double) sample_in_quadrant((double x0, double y0) pos, (double wx, double wy) dir, double lmax, double sigma, Random rng) {
            while (true) {
                var l = sample_length(lmax, rng);
                (double x, double y) = move_point(pos, dir, l);

                var dort = (dir.wy, -dir.wx); // orthogonal to the line direction

                var s = gaussian((0.0, sigma), rng); // could be laplace instead of gaussian

                (x, y) = move_point((x, y), dort, s);
                if (x >= -1.0 && x <= 0.0 && y >= 0.0 && y <= 1.0) // acceptance/rejection
                    return (x, y);
            }
        }

        public static (double, double) sample_in_plane((double x, double y) pos, (double wx, double wy) dir, double lmax, double sigma, Random rng) {
            (double x, double y) = sample_in_quadrant(pos, dir, lmax, sigma, rng);

            if (rng.NextDouble() < 0.25)
                return (x, y);

            if (rng.NextDouble() < 0.5) // reflection over X
                return (x, -y);

            if (rng.NextDouble() < 0.75) // reflection over Y
                return (-x, y);

            return (-x, -y); // reflection over X&Y
        }

        static void Main(string[] args) {
            var rng = new Random(32345);

            var L = 0.5 * SQRT_5 + 0.5 / SQRT_5; // sampling length, BIGGER THAN JUST A SEGMENT IN THE QUADRANT
            (double x0, double y0) pos = (-1.0, 0.0); // initial position
            (double wx, double wy) dir = (2.0 / SQRT_5, 1.0 / SQRT_5); // directional cosines, wx*wx + wy*wy = 1
            double sigma = 0.2; // that's a value to play with

            // last rejection stage
            (double x, double y) pt;
            while(true) {
                pt = sample_in_plane(pos, dir, L, sigma, rng);

                if (pt.x < 0.5) // reject points in the red area, accept otherwise
                    break;
            }
            Console.WriteLine(String.Format("{0} {1}", pt.x, pt.y));
        }
    }
}