Why is my implementation of the parking lot test for random number generators producing bad results?

Question

I'm trying to write an implementation of the parking lot test for random number generators. Here are the sources that I'm getting my information about the test from: Intel math library documentation and Page 4 of this paper along with the phi function for probability density listed here.

I wrote an implementation of the test in C#. It uses a 100x100 grid whose values are initially set to null. I then use the random number generator to generate random integers for x and y. If that index of the grid and it's neighbors are empty, that index gets set to 1. Otherwise, nothing happens because there was a "crash".

I ran it using C# System.Random generator. I don't believe the results are correct because I always get very near 3079 points parked, which is about 500 short of the average I'm supposed to get. It's also yields a p-value of 2.21829146215425E-90.

My code is below. Does anyone have any experience with this or can anyone see something that I might be doing incorrectly in my implementation? Any help would be greatly appreciated.

  private void RunParkingLotTest()
    {
        points = new int?[100,100];
        int parked = 0;

        for (int i = 0; i < 12000; i++)
        {
            int x = random.Next(100);
            int y = random.Next(100);

            if (IsSafeToPark(x, y))
            {
                points[x, y] = 1;
                parked++;
            }

        }
        Console.WriteLine("Parked: " + parked + "\nP value: " + PhiFunction((parked-3523)/21.9));
    }

    private bool IsSafeToPark(int x, int y)
    {
        return PointIsEmpty(x, y) 
            && LeftOfPointIsEmpty(x, y) 
            && RightOfPointIsEmpty(x, y) 
            && BelowPointIsEmpty(x, y) 
            && AbovePointIsEmpty(x, y);
    }

    private bool AbovePointIsEmpty(int x, int y)
    {
        if (y == 99)
        {
            return true;
        }
        else
            return points[x, y + 1] == null;
    }

    private bool BelowPointIsEmpty(int x, int y)
    {
        if (y == 0)
        {
            return true;
        }
        else
            return points[x, y - 1] == null;
    }

    private bool RightOfPointIsEmpty(int x, int y)
    {
        if (x == 99)
        {
            return true;
        }
        else
            return points[x + 1, y] == null;
    }

    private bool LeftOfPointIsEmpty(int x, int y)
    {
        if (x == 0)
        {
            return true;
        }
        else
            return points[x - 1, y] == null;
    }

    private bool PointIsEmpty(int x, int y)
    {
        return points[x, y] == null;
    }

    private double PhiFunction(double x)
    {
        //ϕ(x) = (2π)−½e−x2/2

        return ((1 / Math.Sqrt(2 * Math.PI)) * Math.Exp(-(Math.Pow(x, 2)) / 2));
    }

edit - The problems with my original implementation were

• I was plotting squares instead of disks
• I only plotted points at integer values. I should have used decimal values instead.
• As a result of the above two, I needed to change my distance check

Thanks to Chris Sinclair and mine z for help in figuring this out. The final code is posted below.

Answer 1

I'm going to take a stab at this, and admittedly, I have not attempted any such test, so forgive me if I'm way off. In general though, the .NET Random implementation is pretty good and I've never had issues with it, so I wouldn't suspect that initially especially since you're properly reusing the same instance instead of creating new ones.

Reading from the parking.pdf, and from the Intel documentation, it seems that they're using discs, and compute the distance from their centre points. Your implementation is using squares (array of 1 distance between spots) and thus ignoring diagonals.

From the pdf:

If disks were being used, the distance between the particles r = p(x(i) − z)2 + (y(i) − z)2 would need to be less than or equal to one. Does it matter whether one uses disks or squares? An indication of the importance of which geometric figure is parked can be obtained by comparing the area occupied by a square of side 1.0 to the area of a disk of diameter 1.0. The ratio of the areas, disk to square, is π/4. Therefore, it would be anticipated that more disks could be placed in a box than squares in the same number of tries.

And the Intel doc:

The test assumes a next random point (x, y) successfully ”parked”, if it is far enough from every previous successfully ”parked” point. The sufficient distance between the points (x1, y1) and (x2, y2) is min(|x1 - x2|,|y1 - y2|) > 1.

I'm guessing that the π/4 disk to square ratio and the differences between how many discs can fit vs squares might be why you're seeing a different number. (although right now I fail to see a directly relationship between 3523 and 3070 and π/4. 3523 * π/4 = 2767, which is close, but I'm sure if there's a relationship it's slightly more complex than just simple multiplication.)

Not a great answer, but my best guess.

EDIT: Interestingly enough, I did a quick implementation using discs with 1 unit diameter and getting results around 4000 parked. So maybe there's a bit more to this than my untrained self can grasp (or maybe .NET's Random doesn't pass the test?) Anyway, here's my disc implementation:

List<Point> parkedCars = new List<Point>();
Random random = new Random();

void Main()
{
    int parked = 0;

    for (int i = 0; i < 12000; i++)
    {
        double x = random.NextDouble() * 100;
        double y = random.NextDouble() * 100;

        Point pointToPark = new Point(x, y);

        if (IsSafeToPark(pointToPark))
        {
            parkedCars.Add(pointToPark);
            parked++;
        }

    }
    Console.WriteLine("Parked: " + parked);
}

private bool IsSafeToPark(Point pointToPark)
{
    //make sure it's "inside" the box
    if (pointToPark.X < 0.5 || pointToPark.X > 99.5
        || pointToPark.Y < 0.5 || pointToPark.Y > 99.5)
        return false;

    if (parkedCars.Any(p => Distance(pointToPark, p) <= 1))
        return false;

    return true;
}

private double Distance(Point p1, Point p2)
{
    return Math.Sqrt((p1.X - p2.X) * (p1.X - p2.X) + (p1.Y - p2.Y) * (p1.Y - p2.Y));
}

Using my likely too simple application of the π/4 ratio yields about 3142. A bit closer, but it seems very incorrect.

EDIT: As @mike z pointed out, my test using directly distance is incorrect. According to the parameters of the test, which I forgot about, just checks that the X and Y distance are greater than 1. Changing my Distance check to:

Math.Max(Math.Abs(p1.X - p2.X), Math.Abs(p1.Y - p2.Y))

Yields a much closer result around 3450, which is pretty close. If I take out my "//make sure it's "inside" the box" check, averaged over 10 tries gets 3531!

So my final, "working" code is:

public struct Point
{
    public double X,Y;

    public Point(double x, double y)
    {
        this.X = x;
        this.Y = y;
    }
}

List<Point> parkedCars = new List<Point>();
Random random = new Random();

void Main()
{
    int parked = 0;

    for (int i = 0; i < 12000; i++)
    {
        double x = random.NextDouble() * 100;
        double y = random.NextDouble() * 100;

        Point pointToPark = new Point(x, y);

        if (IsSafeToPark(pointToPark))
        {
            parkedCars.Add(pointToPark);
            parked++;
        }

    }

    Console.WriteLine("Parked: " + parked);
}

private bool IsSafeToPark(Point pointToPark)
{
    if (parkedCars.Any(p => Distance(pointToPark, p) <= 1))
        return false;

    return true;
}

private double Distance(Point p1, Point p2)
{
    return Math.Max(Math.Abs(p1.X - p2.X), Math.Abs(p1.Y - p2.Y));
}

EDIT: I ran the test 100 times twice, and averaged the results to 3521.29 and 3526.74 respectively. Not sure if this means there still something slightly more to this, but perhaps this is just indicative of rounding or floating point precision differences between .NET and Fortran.