Drunkard's Walk

[dutch]

I thought the average distance the drunk would walk is the sqrt of the number of steps. However with 10000 steps I get a distance of about 90, whereas sqrt(10000) is 100.

The answer I get also depends on the number of dimensions and whether the steps are orthogonal or not. Actually 2d_orth and 2d_all are pretty much the same, so maybe that doesn't matter. Hopefully someone can complete the 3d_all mode.

1d:      80
2d_orth: 88
3d_orth: 92
2d_all:  89
3d_all:  ??

#include <iostream>
#include <random>
#include <cstdlib> // abs, atoi
#include <cmath>   // sqrt
using namespace std;

const double Pi = 3.14159265359;
default_random_engine Rnd(random_device{}());

void label(const char *s, int reps, int steps) {
    cout << s << " (reps " << reps << ", steps " << steps << "): ";
}

// 1d: Gives an answer of about 80.
void mode_1d(int reps, int steps) {
    label("1d", reps, steps);
    uniform_int_distribution<int> direction(0, 1);
    int sum = 0;
    for (int rep = 0; rep < reps; ++rep) {
        int n = 0;
        for (int i = 0; i < steps; ++i)
            n += (direction(Rnd) ? 1 : -1);
        sum += abs(n);
    }
    cout << (double)sum / reps << '\n';
}

// 2d_orthogonal: Gives an answer of about 88.
void mode_2d_orthogonal(int reps, int steps) {
    label("2d_orth", reps, steps);
    uniform_int_distribution<int> direction(0, 1);
    double sum = 0;
    for (int rep = 0; rep < reps; ++rep) {
        int x = 0, y = 0;
        for (int i = 0; i < steps; ++i)
            if (direction(Rnd)) x += (direction(Rnd) ? 1 : -1);
            else                y += (direction(Rnd) ? 1 : -1);
        sum += sqrt(x * x + y * y);
    }
    cout << sum / reps << '\n';
}

// 3d_orthogonal: Gives an answer of about 92.
void mode_3d_orthogonal(int reps, int steps) {
    label("3d_orth", reps, steps);
    uniform_int_distribution<int> direction(0, 1), axis(0, 2);
    double sum = 0;
    for (int rep = 0; rep < reps; ++rep) {
        int x = 0, y = 0, z = 0;
        for (int i = 0; i < steps; ++i)
            switch (axis(Rnd))
            {
            case 0: x += (direction(Rnd) ? 1 : -1); break;
            case 1: y += (direction(Rnd) ? 1 : -1); break;
            case 2: z += (direction(Rnd) ? 1 : -1); break;
            }
        sum += sqrt(x * x + y * y + z * z);
    }
    cout << sum / reps << '\n';
}

// 2d_all_directions: Gives an answer of about 89
void mode_2d_all_directions(int reps, int steps) {
    label("2d_all", reps, steps);
    uniform_real_distribution<> direction(0, 2 * Pi);
    double sum = 0;
    for (int rep = 0; rep < reps; ++rep) {
        double x = 0, y = 0;
        for (int i = 0; i < steps; ++i) {
            double angle = direction(Rnd);
            x += cos(angle);
            y += sin(angle);
        }
        sum += sqrt(x * x + y * y);
    }
    cout << sum / reps << '\n';
}

/*
I'm not sure how to do this one.

// 3d_all_directions: Gives an answer of about 
void mode_3d_all_directions(int reps, int steps) {
    label("3d_all", reps, steps);
    uniform_real_distribution<> direction(0, 2 * Pi);
    double sum = 0;
    for (int rep = 0; rep < reps; ++rep) {
        double x = 0, y = 0, z = 0;
        for (int i = 0; i < steps; ++i) {
            double angle1 = direction(Rnd), angle2 = direction(Rnd);
            x += ???;
            y += ???;
            z += ???;
        }
        sum += sqrt(x * x + y * y + z * z);
    }
    cout << sum / reps << '\n';
}
*/

int main(int argc, char **argv) {
    int mode = 0, reps = 1000, steps = 10000;
    if (argc >= 2) mode  = atoi(argv[1]);
    if (argc >= 3) reps  = atoi(argv[2]);
    if (argc >= 4) steps = atoi(argv[3]);
    switch (mode) {
    case 0: mode_1d               (reps, steps); break;
    case 1: mode_2d_orthogonal    (reps, steps); break;
    case 2: mode_3d_orthogonal    (reps, steps); break;
    case 3: mode_2d_all_directions(reps, steps); break;
    //case 4: mode_3d_all_directions(reps, steps); break;
    }
}

[helios]

Something is wrong. How can mode_1d() return 80 after 10000 steps if at each step you added or subtracted 1 with a fair probability? The result should tend towards zero. It sounds like your generator has a bias of 0.8% in either direction.

[dutch]

I don't believe the tendency would be towards zero. As soon as you move away from 0, that position essentially becomes the new zero. There's no "memory" that would tend it to prefer to move back towards zero instead of away from it. You're thinking about the ratio of 0's and 1's tending towards 0.5, which is true. But the absolute difference tends to increase with more "flips of the coin".

Run the program yourself. You'll get the same result. It's best to run it with 10000 reps (of the default 10000 steps) instead of the default 1000. So run it like ./a.out 0 10000

https://en.wikipedia.org/wiki/Random_walk

[Ganado]

I think the ~80 value is representative of the variance, not the average value. If the first run gave a slightly negative answer in total, and the second run gave a slightly positive answer in total, then the sum of the absolute values of those will non-zero. Increasing the number of steps increases the variance. (Sorry, I know this doesn't answer the actual question, I just wanted to throw my two cents in.)

[helios]

Never mind. I was thinking of the average position, not the average translation. Since the final position has a normal distribution with its peak at zero, both sides of the distribution would tend to cancel out, leaving zero.

[Ganado]

dutch, you want to be using spherical coordinates.
https://en.wikipedia.org/wiki/Spherical_coordinate_system

But if you want the points in the volume of the sphere or on the surface of it to be uniformly distributed, you can't just use the normal spherical coordinates (uniform polar angle and azimuthal angle), you have to manipulate it a bit. (That's what I didn't understand at first in the following thread.)
See this thread: http://www.cplusplus.com/forum/windows/262648/

My post in that thread is wrong. Look at Thomas Huxhorn's and lastchance's replies.
Specifically, the getPoint function.

[dutch]

I was thinking of the average position, not the average translation.

Yeah. Since I'm adding the absolute values, it doesn't tend to zero.

I was wondering how to calculate the expected translation and why it was 80. I vaguely remember that it should be the square root of the number of steps, which would be 100 in this case.

[dutch]

you want to be using spherical coordinates.

Thanks ganado. I see it's trickier than I thought. The wolfram page shows multiple ways.
Using one of them gives an answer of 92, the same as the 3d_orthogonal mode.

Probably not the most efficient way since it throws away a lot of points, but:

// 3d_all_directions: Gives an answer of about 92 for 10000 steps.
void mode_3d_all_directions(int reps, int steps) {
    label("3d_all", reps, steps);
    uniform_real_distribution<> dir(-1, 1);
    double sum = 0;
    for (int rep = 0; rep < reps; ++rep) {
        double x = 0, y = 0, z = 0;
        for (int i = 0; i < steps; ) {
            double a = dir(Rnd), b = dir(Rnd);
            double s = a * a + b * b;
            if (s >= 1.0) continue;
            double t = 2 * sqrt(1 - a * a - b * b);
            x += a * t;
            y += b * t;
            z += 1 - 2 * s;
            ++i;
        }
        sum += sqrt(x * x + y * y + z * z);
    }
    cout << sum / reps << '\n';
}

[againtry]

https://www.feynmanlectures.caltech.edu/I_06.html
The famous Feynman lectures. Fig 6-2 and the random walk section discuss and explain the significance of sqrt(N)

In short, the experiment for any particular number of steps needs to be run and re-run many times to get an ‘average’, and obviously the more steps, the greater the ‘average’ is.

My guess is there wouldn’t be much difference between the results even though the no. of dimensions being considered varies, at least not without explanation.

[dutch]

Thanks againtry. Good to see a Feynman reference.

So the sqrt of the average of the squares of the distances (root mean square, RMS) is sqrt(N). Interesting.

But the average of the absolute values isn't sqrt(N). I was able to calculate the expected value of the distance travelled in 10000 steps as 79.786..., which is approximately what the program was getting. The analysis involved Pascal's triangle to determine the number of times a particular result would come up. Ultimately it led to the following bc script:

n = 10000
x = 2; d = 2;
for (i = 2; i <= n; ++i) {
    x *= 2;
    if (i % 2 == 1) {
        x += x / d;
        d += 2;
    }
}
scale = 9
x / (2 ^ n)

So, although the RMS is a (possibly better) indication of the average distance, it isn't actually the average distance. I wonder why it's used?

Here's the 1d version using the RMS.

// Usage: drunk [ REPS [ STEPS ] ]

#include <iostream>
#include <random>
#include <cmath>
using namespace std;

int main(int argc, char **argv) {
    default_random_engine Rnd(random_device{}());
    int reps = 1000, steps = 10000;
    if (argc >= 2) reps  = atoi(argv[1]);
    if (argc >= 3) steps = atoi(argv[2]);
    cout << "reps " << reps << ", steps " << steps << ": ";
    uniform_int_distribution<int> direction(0, 1);
    long sum = 0;
    for (int rep = 0; rep < reps; ++rep) {
        int n = 0;
        for (int i = 0; i < steps; ++i)
            n += (direction(Rnd) ? 1 : -1);
        sum += (long)n * n;
    }
    cout << sqrt((long double)sum / reps) << '\n';
}

[againtry]

Like the Buffon Needle problem to estimate PI, random walk isn't the best way to get the sqrt of a number. But run enough tests and you get an engineer's 'near enough is good enough' answer. Put the distances in a plot and, I haven't tried it, and you'll get something akin to fig 6-2 which is explained via Pascal and binomials blah blah.

#include <iostream>
#include <cmath>
#include <random>
#include <vector>

struct Point
{
    double x;
    double y;
    
    double square_distance()
    {
        return x*x + y*y;
    }
    
    void move(int direction)
    {
        double delta_x{0};
        double delta_y{0};
        
        switch( direction )
        {
            case 1:
                delta_x = 1;
                break;
                
            case 2:
                delta_y = 1;
                break;
            case 3:
                delta_x = -1;
                break;
            case 4:
                delta_y = -1;
                break;
        }
        x += delta_x;
        y += delta_y;
    }
    
};


int main()
{
    int NO_MOVES{100};
    int NO_TESTS{10000};
    
    std::random_device rd;
    std::mt19937 gen(rd());
    std::uniform_int_distribution<> distrib(1, 4);
    
    Point point;
    
    double total{0};
    
    for(int j = 0; j < NO_TESTS; j++)
    {
        point.x = 0;
        point.y = 0;
        
        for (int i = 0; i < NO_MOVES; ++i)
        {
            point.move( distrib(gen) );
            
        }
        total += point.square_distance();
    }
    std::cout << "RMS distance travelled: " << sqrt(total/NO_TESTS) << '\n';
    
    return 0;
}

RMS distance travelled: 9.90872
Program ended with exit code: 0

[againtry]

And,

int NO_MOVES{10000};
    int NO_TESTS{10000};

gives,

RMS distance travelled: 100.586
Program ended with exit code: 0

Time now to get drunk as a skunk, as they say.