How to use the method of Monte Carlo to search for the limit probabilities

Question

I want to understand how to use the method of Monte Carlo to search for the limit probabilities of the some system S.

For example:

    S0  S1  S2  S3
S0  0.1 0.9 0   0
S1  0   0.2 0.3 0.5
S2  0.2 0.1 0.5 0.2
S3  0.5 0   0.4 0.1

As I understand the method we need to generate a some number (x) and then compare the probability:

if x 
   0   <= x < 0.1 => S0 -> S0
   0.1 <= x < 0.9 => S0 -> S1
   0.9 <= x < 0.9 => S0 -> S2
   0.9 <= x < 0.9 => S0 -> S3
   0.9 <= x < 1   => S0 -> S4

when S4 - limit (border)

Similarly for other states.

Following this approach I can count the number of transitions:

  static double[] SimpleMonte(double[][] a, int iter = 1)
    {
        var n = a.GetLength(0);

        var p =
            a
            .Select(x => x.Select((_, i) => x.Take(i + 1).Sum()).ToArray())
            .ToArray();

        Random rand = new Random();

        double[] X = new double[n];
        for (int x = 0; x < n; x++)
        {
            double count = 0;
            for (int i = 0; i < iter; i++)
            {
                int row = x;
                bool notG = true;
                Console.Write("{0} -> ", row);
                while (notG)
                {

                    var e = rand.NextDouble();
                    Console.Write("({0})", Math.Round(e, 2));
                    bool ch = false;
                    for (int j = 0; j < n - 1; j++)
                    {
                        if (p[row][j] <= e && e < p[row][j + 1])
                        {
                            row = j + 1;
                            ch = true;
                            break;
                        }
                    }
                    if (!ch)
                        notG = false;
                    else
                    {
                        Console.Write("{0} -> ", row);
                        count++;
                    }
                }
                Console.WriteLine();
            }
            X[x] = count / iter;
        }
        return X;
    }

https://dotnetfiddle.net/nJF5sm

I will be glad to hear a hint on what to do to solve this problem.

Can you clarify your problem? Do you want to solve `S x = x`, i.e. find the stationary setting? — Nico Schertler, Nov 27 '16 at 15:41

score 2 · Accepted Answer · edited Nov 27 '16 at 16:58

In a practical sense, the best way to find a limit of the a system like that is repeated squaring of the matrix until the entries converge. This works because it is a stochastic matrix (sum of each row equals 1). When I tried it I got an answer of:

S0        S1        S2        S3
0.1939252 0.2593458 0.3294393 0.2172897

which gives the average probabilities that you will be in a particular state.

To use a Monte Carlo method, you should generate the random numbers like you have done and keep count of the transitions. Then the average probability that you're in a certain state is

(Amount of Transitions to State S)/(Total Transitions)

as your Total Transitions gets large enough.

In the code you provided, if you keep increasing the size of your iter variable (and it does need to be relatively large) the last four lines of your output should converge to the numbers above. I hope that helps.

In the original code there was a mistake that prevented a transition to the initial state. That is the correct version:

static double[] SimpleMonte(double[][] a, int iter = 10000)
    {
        var n = a.GetLength(0);

        var p =
            a
            .Select(x => x.Select((_, i) => x.Take(i + 1).Sum()).ToArray())
            .ToArray();

        Random rand = new Random();
        double[] X = new double[n];
        int row = rand.Next(n);
        for (int i = 0; i < iter; i++)
        {
            var e = rand.NextDouble();
            X[row]++;
            if (e < p[row][0])
                row = 0;
            else
                for (int j = 0; j < n - 1; j++)
                {
                    if (p[row][j] <= e && e < p[row][j + 1])
                    {
                        row = j + 1;
                        break;
                    }
                }

        }
        return X.Select(x => x / iter).ToArray();
    }

you know, I took a break and came to this conclusion myself. Thank you :) — Ev_Hyper, Nov 27 '16 at 16:23

How to use the method of Monte Carlo to search for the limit probabilities

1 Answers1